Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Towards ultrasound-based intraoperative dosimetry for prostate brachytherapy Wen, Xu 2010

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

Item Metadata

Download

Media
24-ubc_2010_fall_wen_xu.pdf [ 31.15MB ]
Metadata
JSON: 24-1.0071078.json
JSON-LD: 24-1.0071078-ld.json
RDF/XML (Pretty): 24-1.0071078-rdf.xml
RDF/JSON: 24-1.0071078-rdf.json
Turtle: 24-1.0071078-turtle.txt
N-Triples: 24-1.0071078-rdf-ntriples.txt
Original Record: 24-1.0071078-source.json
Full Text
24-1.0071078-fulltext.txt
Citation
24-1.0071078.ris

Full Text

Towards Ultrasound-based Intraoperative Dosimetry for Prostate Brachytherapy by Xu Wen  B.Eng., Tsinghua University, China, 1999 M.A.Sc., Tsinghua University, China, 2002  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Electrical and Computer Engineering)  The University of British Columbia (Vancouver) July, 2010 c Xu Wen, 2010 ⃝  Abstract Prostate brachytherapy is a widely used treatment of localized prostate cancer. Intra-operative dose feedback would bring many benefits to patients and healthcare practitioners. Detection of brachytherapy seeds and segmentation of prostate boundaries play key roles in dosimetry for prostate brachytherapy. However, seed detection and prostate segmentation using conventional B-mode transrectal ultrasound still remains a challenge for prostate brachytherapy, mainly due to the small size of brachytherapy seeds in the relatively low-quality B-mode ultrasound images and due to the poor contrast between the prostate gland and surrounding tissues, speckle noise, shadowing and refraction artifacts. In this thesis, a new method called the reflected power imaging is presented to enhance the visibility and imaging of implanted seeds. It directly measures the reflected energy of ultrasound radio-frequency signals without logarithmic compression. Based on this method, we propose a new solution for brachytherapy seed detection in a 3D reflected power image computed from ultrasound radio-frequency signals, instead of conventional B-mode images. Then implanted seeds are segmented in 3D local search spaces that are determined by a priori knowledge, e.g. needle entry points and seed placements in a pre-operative dosimetry plan. Needle insertion tracks are also detected locally by using the Hough Transform. Experimental results show that the proposed solution works well for seed localization in a tissue-equivalent ultrasound prostate phantom implanted according to a realistic treatment plan with 136 seeds from 26 needles. As the prostate is a firm organ relative to surrounding tissues, elastography is a potential imaging modality for the guidance of prostate brachytherapy. A dynamic ultrasound elastography method named vibro-elastography can provide more complete dynamic tissue description in terms of transfer functions and coherence functions. In this thesis, we develop fast computational algorithms and programs to implement vibro-elastography imaging in real time. Phantom experiments demonstrate that the vibro-elastography techniques produce stable and operatorindependent strain images with high contrast-to-noise ratio. Furthermore, the software for a 3D vibro-elastography imaging system has been designed, implemented and used in the data collection. Over 15 patients have been scanned at the British Columbia Cancer Agency, Vancouver Centre, and the results are encouraging.  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  List of Tables  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  vi  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  vii  List of Abbreviations  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xii  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xiii  1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  1  Abstract  Acknowledgments  1.1  1.2  Brachytherapy Seed Detection  . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  4  1.1.1  Fluoroscopic Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  4  1.1.2  Ultrasound Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  5  Ultrasound Elastography  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  1.2.1  Tissue Motion Estimation  1.2.2  Vibro-Elastography  1.2.3  8  . . . . . . . . . . . . . . . . . . . . . . . . . . .  9  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  10  Clinical Applications of Elastography . . . . . . . . . . . . . . . . . . . . .  12  1.3  Prostate Boundary Segmentation  . . . . . . . . . . . . . . . . . . . . . . . . . . .  12  1.4  Research Goals and Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  15  1.5  Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  15  2 Enhancement of Seed Visibility in 2D Ultrasound . . . . . . . . . . . . . . . . .  17  2.1  2.2  Candidate Methods  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  17  2.1.1  Reflected Power of RF Signals . . . . . . . . . . . . . . . . . . . . . . . . .  17  2.1.2  Squared Envelope of RF Signals . . . . . . . . . . . . . . . . . . . . . . . .  20  2.1.3  Decompressed B-mode Images . . . . . . . . . . . . . . . . . . . . . . . . .  21  Image Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  21  2.2.1  Sector Reformatting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  21  2.2.2  Image Enhancement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  22  2.2.3  Beam Steering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  22  2.3  Experiments  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  23  2.4  Results  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  25  2.5  Discussion  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  27 iii  3 Brachytherapy Seed Detection Using 3D Ultrasound 3.1  3D TRUS Imaging System  3.2  3D Reflected Power Imaging  3.3  3.4  . . . . . . . . . . . . . . .  35  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  35  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  36  3.2.1  2D Reflected Power Imaging  . . . . . . . . . . . . . . . . . . . . . . . . . .  36  3.2.2  3D Image Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . .  36  Seed Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  37  3.3.1  Volume Cropping  37  3.3.2  Seed Voxel Grouping  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  39  3.3.3  Seed Cluster Squeezing . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  39  3.3.4  Non-seed Artifact Removal . . . . . . . . . . . . . . . . . . . . . . . . . . .  39  3.3.5  Needle Insertion Track Detection  . . . . . . . . . . . . . . . . . . . . . . .  39  3.3.6  Seed Matching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  41  3.3.7  Result Checking against the Pre-Plan . . . . . . . . . . . . . . . . . . . . .  41  3.3.8  False Positive Removal and/or Missing Seed Identification  . . . . . . . . .  42  3.3.9  Seed Localization  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  42  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  Fluoroscopic Tomosynthesis  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  42  3.4.1  Seed Segmentation in 2D Fluoroscopic Images  . . . . . . . . . . . . . . . .  43  3.4.2  Back Projection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  43  3.4.3  Phantom Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  45  3.5  Experiments  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  47  3.6  Results  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  49  3.7  Discussion  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  57  4 Fast Computation Algorithms and Programs for Vibro-Elastography . . . . .  62  4.1  Transfer Function Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  62  4.1.1  Moving Windows  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  63  4.1.2  Accumulative Windows . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  64  4.1.3  Average of Transfer Functions  . . . . . . . . . . . . . . . . . . . . . . . . .  65  4.2  Power Strain Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  66  4.3  Application Programs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  66  4.3.1  RF Data Analysis Program . . . . . . . . . . . . . . . . . . . . . . . . . . .  66  4.3.2  Real-time Data Acquisition Program  . . . . . . . . . . . . . . . . . . . . .  69  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  71  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  72  4.4  Experiments  4.5  Results  4.6  Discussion  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  73  5 Vibro-Elastography Imaging of the Prostate Region . . . . . . . . . . . . . . . .  77  5.1  5.2  VE System Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  77  5.1.1  Ultrasound Machine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  77  5.1.2  Actuation System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  77  3D VE Imaging 5.2.1  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  2D Signal/Image Processing  . . . . . . . . . . . . . . . . . . . . . . . . . .  79 79 iv  5.2.2  3D VE Image Reconstruction  . . . . . . . . . . . . . . . . . . . . . . . . .  81  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  81  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  82  5.3  Experiments  5.4  Results  5.5  Discussion  6 Conclusions  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  84  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  99  6.1  Contributions  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  99  6.2  Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103  Appendices A Ethics Certificate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 B Consent Form for Vibro-Elastography Volume Study . . . . . . . . . . . . . . . 115 C Consent Form for Vibro-Elastography Study in the OR  . . . . . . . . . . . . . 120  D Data Acquisition Protocols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 D.1 Data Acquisition Protocol for Vibro-Elastography Volume Study . . . . . . . . . . 127 D.1.1 Purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 D.1.2 Study location and duration  . . . . . . . . . . . . . . . . . . . . . . . . . . 127  D.1.3 Summary of the procedure and equipment  . . . . . . . . . . . . . . . . . . 127  D.1.4 Suggested equipment setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 D.1.5 Data Acquisition Steps  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128  D.2 Data Acquisition Protocol for Vibro-Elastography Study in the OR D.2.1 Study location and duration  . . . . . . . . 129  . . . . . . . . . . . . . . . . . . . . . . . . . . 129  D.2.2 Equipment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 D.2.3 Suggested equipment setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 D.2.4 Data Acquisition Steps  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129  E Communication Protocol between the Data Acquisition Program and VE Actuation System  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131  E.1 Writing Shared Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 E.2 Reading Shared Files  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131  F Recipes for Making Phantoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 F.1 Equipment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 F.2 Making agar phantoms  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133  F.3 Making gelatin phantoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133  v  List of Tables 2.1  Contrast and CNR between all the seed regions and the background in the B-mode, decompressed B-mode, reflected power and squared envelope images, in the three transverse scan planes in Figs. 2.5 (a)-(d). . . . . . . . . . . . . . . . . . . . . . . .  2.2  Statistics of the absolute errors between the seed centroids in the transverse reflected power images and the corresponding CT images. . . . . . . . . . . . . . . . . . . .  2.3  26 26  Running time for computing a frame of decompressed B-mode, reflected power or squared envelope image, respectively, relative to the data acquisition frame period, using the Ultrasonix computer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  3.1  26  Statistics of the seed localization errors by the tomosynthesis reconstruction. The coordinate frame is defined as the same in Fig. 3.1 (b). . . . . . . . . . . . . . . . .  45  3.2  Statistics of the seed localization errors between the 3D ultrasound and CT scans.  51  3.3  Statistics of the seed localization errors between the 3D ultrasound and fluoroscopic tomosynthesis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  3.4  Statistics of the dose curve errors between the 3D ultrasound and CT scan, and between the 3D ultrasound and fluoroscopic tomosynthesis. . . . . . . . . . . . . .  4.1  57  Contrast and CNR between the hard inclusion and the background in the VE and power strain images (Fig. 4.8 (b) and (c)) of the hard inclusion phantom. . . . . .  4.2  57  72  Computational time per frame for power strain imaging and VE using accumulative windows in the phantom experiments. The total time includes the TDPE, PSD and/or TF estimation, and LSE of strain. . . . . . . . . . . . . . . . . . . . . . . .  5.1 5.2  72  Summary of the data acquisition (DAQ) protocols used for VE imaging of the prostate region during the volume study and in the OR prior to the implant. . . .  82  Parameters used for VE image computation offline. . . . . . . . . . . . . . . . . . .  82  vi  List of Figures 1.1  Sketch map of the anatomy of the prostate region and of prostate brachytherapy (reprinted from [107], and no copyright is associated with this image). . . . . . . .  1.2  1  Sketch map of the geometry of TRUS-guided prostate brachytherapy (reprinted from [115]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  3  1.3  Schematic diagram of a motorized fluoroscopic C-arm. . . . . . . . . . . . . . . . .  5  1.4  Schematic diagram of the procedure of taking fluoroscopic images at different rota-  1.5 1.6  tion angles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  6  Schematic diagram of ultrasound elastography or strain imaging. . . . . . . . . . .  9  The linear dynamic tissue model in vibro-elastography, denoted by  Hij ,  and xi and  xj are the time sequences of the displacements of tissue elements i and j at different spatial locations, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7  11  Example of the prostate B-mode TRUS image in the transversal scan plane, where the TRUS probe is on the bottom and is normal to the image. There are three implanted seeds close to the anterior boundary of the prostate in the image. . . . .  2.1  13  Example of dividing an RF scan line into a series of N short segments with equal length and overlap. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  18  2.2  Flow diagram of computing the reflected power in the frequency domain.  . . . . .  19  2.3  Flow diagram of computing the reflected power in the time domain. . . . . . . . .  20  2.4  Prostate phantom made of gelatin and used for the experiments. . . . . . . . . . .  24  2.5  (a) Original B-mode images, (b) reflected power images, (c) squared envelope images, and (d) decompressed B-mode images of the seeds in the phantom, in three different transversal scan planes (from the left to right), where the TRUS probe is on the bottom and normal to the scan planes. . . . . . . . . . . . . . . . . . . . . .  2.6  30  (a) Original B-mode images, (b) reflected power images, (c) squared envelope images, and (d) decompressed B-mode images of the seeds in the phantom, in one sagittal scan planes, where the TRUS probe is on the bottom and parallel to the scan plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  2.7  31  Normalized curves extracted from (a) the original B-mode TRUS image, (b) the corresponding reflected power image, (c) squared envelope image, and (d) decom-  2.8  pressed B-mode image, in a transverse plane (the middle column of Fig. 2.5). . . .  32  CT scan slices of the seeds in the gelatin phantom, in the three transverse planes. .  32  vii  2.9  3-D visualization results of the 24 implanted seeds in the phantom from two different view angles, using Matlab functions, based on a series of sagittal reflected power images acquired when the TRUS probe was slowly rotated around its long axis by hand.  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  33  2.10 Reflected power images of a regular seed (top row) and an echo seed (bottom row) on the surface of the acoustically transparent phantom, at the beam angles from 0◦ to 25◦ away from the perpendicular in 5◦ increments, in the sagittal scan plane. . .  33  2.11 Reflected power images of a regular seed (top row) and an echo seed (bottom row) on the surface of the acoustically transparent phantom, at the beam angles from 5◦ to 10◦ away from the perpendicular in 1◦ increments, in the sagittal scan plane. . .  33  2.12 Reflected power curves of individual seed types as determined by the sum of intensities (SOI) from the sagittal reflected power images, (a) at beam angles from 0◦ to 25◦ in 5◦ increments and (b) at beam angles from 0◦ to 10◦ in 1◦ increments away from the perpendicular. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.13 Reflected power images of a regular seed at the beam angles of (a)  10◦  and (b)  34  15◦ ,  and reflected power images of an echo seed at the beam angles of (c) 10◦ and (d) 15◦ , after performing beam steering, in the sagittal scan plane. . . . . . . . . . . . 3.1  34  (a) Photograph of the 3D TRUS imaging system, and (b) schematic diagram showing the coordinates used and the rotation motor. The motor and cam mounted on the transducer cradle are used to impart radial motion to the TRUS for elastography images. This feature of the system is not used in this chapter. . . . . . . . . . . . .  36  3.2  Flow diagram of seed detection steps. . . . . . . . . . . . . . . . . . . . . . . . . .  38  3.3  Schematic diagram of the representation of lines in 3D spaces. . . . . . . . . . . . .  40  3.4  Schematic diagram of the geometry of the GE OEC 9800 C-arm, where the definitions of (x, y, z), (u, v) and d are the same as in [138]. . . . . . . . . . . . . . . . .  3.5  (a) Original fluoroscopic image of the seed phantom at  0◦ ,  44  and (b) the seed segmen-  tation result in a cropped region of interest. . . . . . . . . . . . . . . . . . . . . . .  45  3.6  Actual seed distribution in the precisely fabricated seed phantom. . . . . . . . . . .  46  3.7  Detected seed distribution in the precisely fabricated seed phantom by the fluoroscopic tomosynthesis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  3.8  CIRS Model 053 tissue-equivalent ultrasound prostate phantom with the implanted seeds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  3.9  46 47  Distribution of the needle entries in the preoperative dosimetry plan for the CIRS phantom. Needle indices are marked on the top-left corner of every needle entry. The number of seeds by each needle is marked inside the contour of each needle entry. The entry point of needle 1 was moved to c4.5 during the implant procedure.  48  3.10 Distribution of the needle entries in the preoperative dosimetry plan for a patient. Needle indices are marked on the top-left corner of every needle entry. The number of seeds by each needle is marked inside the contour of each needle entry. . . . . .  49  3.11 Seed detection result from the 3D CT scan. . . . . . . . . . . . . . . . . . . . . . .  50  3.12 Seed detection result from the fluoroscopic tomosynthesis. . . . . . . . . . . . . . .  50 viii  3.13 (a) Original B-mode TRUS image and (b) reflected power image, of the seeds in the CIRS phantom in one sagittal scan plane, where the TRUS probe is on the bottom and parallel to the scan plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  51  3.14 Detection result of the implanted seeds in the 3D reflected power image, after thresholding. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  52  3.15 Detection result of the implanted seeds in the 3D reflected power image, after seed voxel grouping. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  52  3.16 Detection result of the implanted seeds in the 3D reflected power image, after seed cluster squeezing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  53  3.17 Detection result of the implanted seeds in the 3D reflected power image, after nonseed artifact removal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  53  3.18 Detection result of needle insertion tracks (red lines) by the Hough transform in 3D spaces. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  54  3.19 Detection result of the implanted seeds in the 3D reflected power image, after false positive removal and missing seed identification. . . . . . . . . . . . . . . . . . . . .  54  3.20 Dose curves (at the prescription dose 144.0 Gy) computed from the ultrasound, CT and fluoroscopic seed detection results, from (a) z = 5 mm to (f) z = 30 mm, at intervals of 5 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  55  3.21 Dose curves (at the prescription dose 144.0 Gy) computed from the ultrasound, CT and fluoroscopic seed detection results, from (a) z = 35 mm to (d) z = 50 mm, at intervals of 5 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  56  3.22 Detection result of the implanted seeds in the 3D reflected power image of the patient. 57 3.23 Detection result of the needle insertion tracks (red lines) by the Hough transform in the 3D reflected power image of the patient. . . . . . . . . . . . . . . . . . . . .  58  4.1  Example of overlapping windows used for vibro-elastography. . . . . . . . . . . . .  63  4.2  Flow diagram for the TF/CF estimation with moving windows. . . . . . . . . . . .  64  4.3  Flow diagram for the TF/CF estimation with accumulative windows. . . . . . . . .  65  4.4  Main GUI of the RF data analysis program, which consists of the main frame with the menu and tool bar, the control panel on the left side, and the image client window on the right side. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  67  4.5  Advanced dialog of the RF image properties with the RF data analysis program. .  68  4.6  Main GUI of the real-time data acquisition (DAQ) program, which consists of the main frame with the menu and tool bar, the control panel on the left side, and the image client window on the right side. . . . . . . . . . . . . . . . . . . . . . . . . .  70  4.7  A biplane TRUS transducer mounted on a motorized actuator. . . . . . . . . . . .  71  4.8  (a) B-mode ultrasound image, (b) VE (TF strain) image and (c) power strain image of the hard inclusion phantom, (d) center vertical A-lines extracted from (b) and (c), and (e) PSD curves at different depths along the center vertical A-line of (c). Color coding of the power strain and VE images is HOT as defined in Matlab (please see the index bar). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  73 ix  4.9  (a) B-mode image and (b) VE image of the prostate phantom in a transversal scan plan, (c) VE image of the prostate phantom in a sagittal scan plane and (d) relative compliance along the marked VE A-line showing a change of a factor of 2, (e) TFs and (f) corresponding CFs along the marked A-line are displayed. The motion applied is 4 mm peak-to-peak, band-pass filtered (0.5-10 Hz) white noise. Transfer functions are averaged over 10 seconds, and are computed with reference to the motion in the focal area of the firing crystal at a depth of 25 mm. Color coding is HOT as defined in Matlab. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  5.1  Photo of the Sonix RP PC-based ultrasound machine (Ultrasonix Medical Corp. Burnaby, BC, Canada). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  5.2  74  78  Schematic drawing of the first-generation VE actuation system. The endorectal TRUS probe is mounted in a housing and is vibrated along its axial plane by two motor-driven four-bar linkages, along a direction set by the rotational stepper, with its depth adjusted by the translational stepper. . . . . . . . . . . . . . . . . . . . .  5.3  79  (a) Schematic drawing and (b) photo of the second-generation actuation system. The motor and cam mounted on the transducer cradle are used to impart radial vibration to the TRUS probe for VE imaging in the sagittal scan planes. . . . . . .  5.4  80  (a) - (c) B-mode TRUS images and (d) - (f) VE images of the prostate region of patient 3 in three sagittal planes at 0◦ , 10◦ and 20◦ , (g) the relative compliance curve, (h) four TF curves and (i) the corresponding CF curves along the A-line in (c), at depths of 1, 2, 3 and 4 cm. Color coding is HOT as defined in Matlab (please see the index bar). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  5.5  83  (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 4, in the sagittal planes at -30◦ , -15◦ , 0◦ , 15◦ and 30◦ , from the top down. Color coding is JET as defined in Matlab (please see the index bar). . . . . . . . . . . . . . . . . . . . . . . . . . . .  5.6  86  (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 5, in the sagittal planes at -20◦ , -10◦ , 0◦ , 10◦ and 20◦ , from the top down. . . . . . . . . . . . . . . . . . . . .  5.7  87  (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 6, in the sagittal planes at -10◦ , -5◦ , 0◦ , 5◦ and 10◦ , from the top down. The white arrow in (j) points to the urethra. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  5.8  88  (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 7, in the sagittal planes at -20◦ , -10◦ , 0◦ , 10◦ and 20◦ , from the top down. . . . . . . . . . . . . . . . . . . . .  5.9  89  (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 8, in the sagittal planes at -20◦ , -10◦ , 0◦ , 10◦ and 20◦ , from the top down. . . . . . . . . . . . . . . . . . . . .  90  x  5.10 (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 9, in the sagittal planes at -20◦ , -10◦ , 0◦ , 10◦ and 20◦ , from the top down. The white arrow in (f) points to the calcification. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  91  5.11 (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 10, in the sagittal planes at -20◦ , -10◦ , 0◦ , 10◦ and 20◦ , from the top down. The white arrows in (f) and (n) point to the calcification.  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  92  5.12 (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 11, in the sagittal planes at -10◦ , -5◦ , 0◦ , 5◦ and 10◦ , from the top down. . . . . . . . . . . . . . . . . . . . . . .  93  5.13 (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 12, in the sagittal planes at -10◦ , -5◦ , 0◦ , 5◦ and 10◦ , from the top down. . . . . . . . . . . . . . . . . . . . . . .  94  5.14 (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 13, in the sagittal planes at -20◦ , -10◦ , 0◦ , 10◦ and 20◦ , from the top down. The white arrow in (n) points to the catheter in the urethra. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  95  5.15 (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 14, in the sagittal planes at -30◦ , -15◦ , 0◦ , 15◦ and 30◦ , from the top down. The white arrow in (j) points to the urethra. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  96  5.16 (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 15, in the sagittal planes at -30◦ , -15◦ , 0◦ , 15◦ and 30◦ , from the top down. The white arrow in (j) points to the urethra. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  97  5.17 Average normalized correlation coefficient and average CF curves of patients 4 - 15, excluding the bladder area, at rotation angles between -45◦ and 45◦ in 5◦ increments. 98  xi  List of Abbreviations 2D  Two Dimensional  3D  Three Dimensional  BC  British Columbia  CF  Coherence Function  CNR  Contrast-to-Noise Ratio  CSD  Cross Spectral Density  CT  Computed Tomography  DAQ  Data Acquisition  FFT  Fast Fourier Transform  HIFU  High Intensity Focused Ultrasound  LSE  Least Squares Estimation  MRI  Magnetic Resonance Imaging  MRSI  Magnetic Resonance Spectroscopic Imaging  PSD  Power Spectral Density  RF  Radio Frequency  RMS  Root-Mean-Square  ROI  Region of Interest  RTD  Real-Time Dosimetry  TDE  Time Domain cross-correlation Estimation  TDPE  Time Domain cross-correlation with Prior Estimates  TF  Transfer Function  TRUS  Transrectal Ultrasound  TUUS  Trans-urethral Ultrasound  VE  Vibro-Elastography  VOI  Volume of Interest  xii  Acknowledgments First of all, I would like to thank my supervisors Dr. Tim Salcudean and Dr. Peter Lawrence for their continuing support during this project, for giving novel ideas and directions to my research, and for helping me with insights and methodologies toward a useful solution. I would gratefully appreciate Dr. Tim Salcudean for designing and building the vibro-elastography actuation systems. We could not have collected the number of patient data without the actuation systems. I also want to thank Dr. Rob Rohling for his help and advices on my research project. I would like to acknowledge the support of Dr. James Morris, Dr. Ingrid Spadinger, Dr. Mira Keyes, Dr. Michael McKenzie, Dr. Tom Pickles and Mr. Nick Chng, in the British Columbia Cancer Agency, Vancouver Centre, who assisted me for obtain patients’ consent for participating our study and for patient data acquisition. Dr. Morris was instrumental in motivating this work and was performed seed implantation to a tissue-equivalent ultrasound prostate phantom. Dr. Spadinger and Mr. Chng generated the pre-operative dosimetry plan for the prostate phantom. Mr. Chng was helping me with patient data collection in the operating room. I would also like to thank Mehdi Moradi and Sara Mahdavi for their support and help with patients’ ultrasound data collection in the British Columbia Cancer Agency, Vancouver Centre and with patients’ MRI data collection in the UBC hospital. I would also appreciate Danny French for all his help of contact with the physicians and medical physicists, of collecting patient data, and of making phantoms, in the past five years. My phantom and patient data collection would not have been successful without the assistance of Andrew Thompson, Simon Bachmann and Honza Vohradsky who built and maintained the VE actuation systems, of Reza Zahari-Azar and Neerav Patel who initially developed the programs for real-time ultrasound data acquisition and tissue motion tracking, of Thomas Diego Prananta who helped me with the mechanical and electrical device issues, and of Julio Lobo and Michael Tien who developed the back projection program for fluoroscopic tomosynthesis. Moreover, I also appreciate the help, suggestions, brainstorms and conversations with everyone in the Robotics and lab. Thank you very much for your kindness and contributions, which have made the Robotics and Control Lab a great open place to work. Additionally, I would like to acknowledge the support of this research by the Natural Sciences and Engineering Research Council (NSERC) of Canada and by the National Institutes of Health (NIH) of the United States. Finally, this thesis is dedicated to my wife and to my parents.  xiii  Chapter 1  Introduction The prostate is a compound tubuloalveolar exocrine gland in the male reproductive system, as shown in Fig. 1.1. Prostate cancer is the most frequently diagnosed cancer and one of the leading causes of cancer death in North American men. According to the Canadian Cancer Society, Canadian men would continue to see prostate cancer as the most commonly diagnosed cancer with an estimated 25,500 new cases [1], and prostate cancer is third in mortality, causing 4,400 deaths in 2009 [1]. In the report of the American Cancer Society, an estimated 192,280 new cases of prostate cancer would occur in the United States during 2009 [2], and prostate cancer is the second-leading cause of cancer death with an estimated 27,360 deaths [2]. At present, the only well-established risk factors for prostate cancer are age, race/ethnicity, and family history of the disease. About 63% of all prostate cancer cases are diagnosed in men aged 65 and older [2].  Figure 1.1: Sketch map of the anatomy of the prostate region and of prostate brachytherapy (reprinted from [107], and no copyright is associated with this image).  Early prostate cancer usually does not have symptoms. Patients may experience weak or interrupted urine flow, inability to urinate or difficulty starting or stopping the urine flow, the need to urinate frequently especially at night, blood in the urine, or pain or burning with urination, with more advanced disease [2]. Advanced prostate cancer commonly spreads or metastasizes to 1  the bones, which can cause pain in the hips, spine, ribs, or other areas [2]. However, most of the prostate cancers grow slowly or remain asymptomatic before extensive local growth or metastasis of the tumor happens or until the patients die of other diseases. Prostate cancer is more likely to be curable when the tumor is diagnosed at an early stage (in particular when the tumor is encapsulated within the prostate gland). Although some treatment options may still work in the late stages, early diagnosis and effective treatment of prostate cancer are significant to patients. As the prostate is well protected inside the pelvis (as shown in Fig. 1.1), it is difficult to examine the prostate gland clinically. This work used to be performed by the physical palpation, prostatic fluid inspection, biopsy, and surgery before. At present, commonly-used screening techniques for diagnosis of prostate cancer include the digital rectal examination (DRE), prostate specific antigen (PSA) test, transrectal ultrasound (TRUS), computed tomography (CT), magnetic resonance imaging (MRI), and magnetic resonance spectroscopic imaging (MRSI) [11]. Ultrasound has been the primary imaging modality in prostate-related applications for the following reasons: • It is not inferior to either MRI or CT in terms of diagnostic value [118]; • It can follow anatomical deformations in real time during diagnosis, biopsy and treatment; • It is inexpensive and is easier to use compared with MRI or CT; • It is safe to both patients and healthcare practitioners, because there is no ionizing radiation involved. Various treatment options are available for prostate cancer that is clinically confined to the prostate gland, e.g. radical prostatectomy (surgical excision), electrodessication and curettage (tissue destruction by electric current and removal by scraping with a curette), hormonal therapy, cryotherapy (tissue destruction by freezing), high intensity focused ultrasound (HIFU) therapy, external beam radiation therapy, and permanent trans-perineal brachytherapy. Among these treatment options, permanent prostate brachytherapy, which is also a form of radiation therapy, is a minimally invasive procedure for treating localized prostate cancer. Fig. 1.2 shows a sketch of the geometry of TRUS-guided prostate brachytherapy. In this treatment, a number of small radioactive sources or seeds, typically  125 I  or  103 Pd,  are permanently implanted  into the prostate by brachytherapy needles under the guidance of TRUS, to deliver a sufficient radiation dose to kill the cancerous cells and maintain a tolerable dose to healthy structures, such as the bladder, urethra and rectum [155]. As the control rates of prostate cancer appear to be dose dependent, it is theorized that the higher doses produced by brachytherapy will yield higher control rates than external-beam radiation [35]. Furthermore, the cost and suffering of patients of brachytherapy are also less than external-beam radiation therapy. Thus, brachytherapy may be a superior treatment option of early-stage prostate cancer than surgery and external-beam radiation therapy. Prostate brachytherapy is associated with excellent biochemical control and long-term outcomes for patients with localized prostate tumors [108, 112].  2  Figure 1.2: Sketch map of the geometry of TRUS-guided prostate brachytherapy (reprinted from [115]).  The efficacy of prostate brachytherapy is evaluated by the delivered radiation dose. In order to deliver a high conformal dose safely to the prostate, the radioactive seeds must be positioned accurately within the prostate gland [111, 125, 126]. Clinically relevant delivered radiation dose measures are used for this purpose, such as the D90 (the radiation dose delivered to 90% of the prostate gland, which is predictive of patient outcomes) or V100 (the prostate volume receiving 100% of the prescription dose) [93, 94, 152]. In order to compute delivered dose, the distribution of the brachytherapy seeds with respect to the boundary of the prostate after implantation and knowledge of the radioisotope source type are used. Typically, the prostate boundary is delineated from pre-operative B-mode TRUS images, and the seed distribution is obtained from post-operative CT scans performed within several days or several hours of the implant. Dose errors can be caused by several factors. Published data supports the generally accepted view that prostate segmentation in B-mode TRUS images is not always reliable, whether done manually or automatically. This is especially the case at the apex where the prostate blends into the muscles of the pelvic floor, and at the base where the anterior portion is iso-echoic with the neck of the bladder. Pre-operative to intra-operative registration of TRUS images is done manually with significant human intervention. During each needle insertion, needle forces push the prostate [27, 73] and cause seed placement errors [9, 10]. Prostate swelling and deformation caused by edema, and prostate motion due to repositioning the TRUS transducer (e.g., to avoid the pubic arch) and due to patients’ respiration also generate errors [94, 108]. Thus, there is a need for intra-operative dose feedback to enable a radiation oncologist to intraoperatively adjust the planned seed placement to achieve the appropriate radiation dose and prevent overdosing of healthy tissues. Intra-operative dose computation (real-time dosimetry or RTD) has been advocated in order 3  to ensure better quality control (repeat procedures are sometimes required, especially at new treatment centers) and the ability to modify the plan and treatment intra-operatively [94], when the patient is still under anesthesia. Intra-operative planning with RTD has been shown to achieve better delivered dose in some studies [152]. Although inter-operative planning with RTD increases procedure time, it eliminates the need for the pre-operative and post-operative visits of patients [94, 152]. An overview of the literature in the detection of brachytherapy seeds, ultrasound elastography techniques, and prostate boundary segmentation algorithms is presented in Sections 1.1, 1.2 and 1.3, respectively.  1.1  Brachytherapy Seed Detection  Detection of brachytherapy seeds plays a key role in dosimetry for prostate brachytherapy. Typically, brachytherapy seeds are cylindrical in shape with a diameter of 0.8 mm and length of 4.5 mm. There are three types of seeds currently being used for prostate brachytherapy: • Regular or standard seeds that have a smooth cylindrical surface; • Echo seeds that have a corrugated surface and incorporate the use of alternating ridges and grooves in a regular sinusoidal pattern along the length of the seed except at the two ends [25, 132]; • Linked or stranded seeds in which the regular seeds and spacers are linked together in a stiffened Vicryl suture material [25].  1.1.1  Fluoroscopic Methods  As mentioned above, seed distribution is currently obtained from post-operative CT scans. However, CT is not a practical imaging modality for intraoperative uses. In addition to TRUS, X-ray fluoroscopy is often used in research for prostate brachytherapy. It is the gold standard for visualizing brachytherapy implants without moving the patient. Typically, individual seeds are located in fluoroscopic images that are acquired from three or more perspectives [56, 119, 129, 134, 138, 154]. Figs. 1.3 and 1.4 show the schematic diagrams of a fluoroscopic C-arm and the procedure of taking fluoroscopic images at different rotation angles, respectively. Tubic et al. [133, 134] utilized a simulated annealing algorithm to automatically match seed artifacts in three radiographic films or fluoroscopic images. The matched seed artifacts were backprojected to compute the 3D coordinates of each seed. Todor et al. [129] proposed an operator-free method with projections of seed positions in three isocentric images taken at arbitrary angles. The method also consists of a set of heuristic rules (in a sense, a learning algorithm) that attempts to minimize seed misclassifications [129]. Tutar et al. [138] developed a brachytherapy-specific tomosynthesis reconstruction algorithm with built-in blur compensation. In contrast to other projection-based methods, their method can detect completely overlapping seeds [138]. Singh et al. [119] proposed an optimization algorithm based on geometric and linear programming techniques. Fichtinger et al. [56, 57] formalized seed-matching as a combinatorial optimization problem 4  which has salient features and utilized the Hungarian algorithm for matching and reconstruction of brachytherapy seeds.  Figure 1.3: Schematic diagram of a motorized fluoroscopic C-arm.  Although fluoroscopy is able to provide the spatial distribution of seeds, it cannot register this distribution with respect to the prostate boundary because fluoroscopy does not provide sufficient soft tissue contrast. To compute the delivered dose, a number of methods fusing TRUS and fluoroscopic images to compute dosimetry have been reported [16, 38, 39, 44, 55, 100, 127, 130, 137]. With these methods, prostate contours are delineated in B-mode TRUS images, and brachytherapy seeds are located in fluoroscopic images. However, updating the dosimetry is also time-consuming and inconsistent with the current TRUS-guided brachytherapy procedure because the C-arm must be rotated to several different angles, with the endorectal TRUS probe removed, to acquire the fluoroscopic images. Furthermore, an image registration step involving TRUS-segmented prostate contours is still required in order to locate the determined cloud of seeds relative to the prostate. Costs are increased by the use of X-ray fluoroscopy, and so is exposure to radiation of patients and healthcare providers.  1.1.2  Ultrasound Methods  As TRUS is the predominant imaging modality providing guidance for prostate brachytherapy, an ultrasound-based solution for intra-operative dose computation would offer many benefits. Many attempts have been made for seed localization in B-mode ultrasound images. However, conventional B-mode ultrasound has difficulties in imaging and detection of brachytherapy seeds for the following reasons [34, 48]: 5  Figure 1.4: Schematic diagram of the procedure of taking fluoroscopic images at different rotation angles.  • The size of brachytherapy seeds is quite small relative to the size of other features in tissue; • The brightness of brachytherapy seeds varies depending on its orientation with respect to the incident ultrasound beam; • Visualization of brachytherapy seeds is often hampered by calcifications, speckle, and shadow artifacts inherent to tissue. All the above-mentioned factors result in B-mode TRUS images with low contrast between the implanted seeds and surrounding tissue. Thus, in a study of manual segmentation of seeds in patients, Han et al. [48] concluded that they could not reliably identify a sufficient percentage of the seeds in B-mode TRUS images for accurate dose computation. Xue et al. [148] tried to manually localize linked seeds and spacers with a Hitachi Model EUB6500 high-resolution TRUS probe (Hitachi America, Brisbane, CA). The advantage of linked seeds is that it is not necessary to visualize all the seeds. Seeds that are invisible in B-mode TRUS images can be localized based on their known location in the needle insertion array. Their approach has not been tested on loose seeds, although the approach may be very viable because more and more linked seeds are being used in clinical practice. Davis et al. [25] provided a quantitative measurement of ultrasound backscatter signals from different seed types (regular seeds, echo seeds and stranded seeds) as a function of incident ultrasound beam angles. They found that corrugated seeds produced a more readily identifiable seed image over a large range of seed orientations compared with regular seeds. 6  Yu et al. [149] made use of distributed constant false alarm rate (CFAR) processors and orientation-sensitive morphological filtering. Tissue scatterers were modeled as a multiplicative and signal-dependent K-distribution. The CFAR problem was posed as the detection of a fluctuating target against K-distributed clutter [149]. The orientation-sensitive morphological filters were used to reduce the speckle noise. Thus, the algorithm was mathematically complex and may not work in the clinical environment with interference from multiple seeds. Holmes et al. [51–53] used trans-urethral ultrasound (TUUS) for prostate imaging and seed detection. TUUS can image the prostate from the urethra through the center of the prostate gland, as opposed to TRUS that images the prostate from the posterior boundary. TUUS requires a smaller imaging depth, improves the spatial resolution [51], and reduces the effects of seeds in shadowing of one another, although it is not the standard imaging modality currently used for prostate brachytherapy. Their segmentation framework incorporated a priori information about the size, shape, and orientation of the seeds with fuzzy inference rules [53]. Ding et al. [28, 29] and Wei et al. [145, 146] developed automated seed segmentation methods in 3D B-mode TRUS images. Their approaches needed to perform a subtraction between a prescan before the needle insertion and a post-scan after the needle insertion, and relied on the needle segmentation result in the 3D difference image. Therefore, their methods would be subject to errors in the presence of prostate’s deformation and patient’s motion between the pre-scan and the post-scan. Additionally, the method proposed by Ding et al. [28, 29] also used a tri-bar model and 3D line segment patterns for seed candidate recognition, while the tri-bar model and 3D line segment patterns may be sensitive and affected by seed orientations and speckle noise. In addition to conventional B-mode ultrasound, other ultrasonic techniques have also been tried for seed detection. McAleavey et al. [87–89] modified brachytherapy seeds to be ferromagnetic. The modified seeds were vibrated using an external magnetic field and then imaged using Power Doppler. This approach was tested on modified seeds in a phantom with some success. Nevertheless, ferromagnetic seeds may not be allowed for prostate brachytherapy, and additional magnetic field or equipment is required for vibrating the ferromagnetic seeds. Mitri et al. [90–92] utilized vibro-acoustography to vibrate and detect metal seeds. They tested the method with an agar phantom in a water tank. Two slightly shifted ultrasound beams were focused on each seed which causes an oscillating force. Because the seeds are inhomogeneous inside the phantom, a variation in the acoustic emission field is observed and used to detect the seeds [91, 92]. In the experiments, the best excitation vibration frequencies were determined either by calculating fundamental resonance frequencies for the seeds or by the experimental optimal resonance frequency of the gel. While showing promise, the need for additional instrumentation makes this approach clinically viable in the long term. Feleppa et al. [34] used signature and elastography techniques with ultrasound radio frequency (RF) signals. The signature method performs 2D cross-correlation between the “signature” signals of a seed in an acoustically transparent medium and the signals of the entire scan plane of the seeds in a phantom. The seed signature varies with different beam angles and with angles between the long axis of the seed and the beam. Correlation analysis using RF signals is also sensitive to changes in phase. The elastography technique requires two or more frames of RF data under varying 7  external compression to calculate one result strain image, and it often fails to detect implanted seeds placed behind each other facing the incident beam as the motion tracking algorithms used for elastography may no longer work behind strong reflectors. Mamou et al. [84] proposed an algorithm based on singular spectrum analysis (SSA), which utilizes eigenvalues derived from the diagonalized correlation matrix of envelope-detected RF echo signals to yield a P value indicative of the likelihood of a seed-specific repetitive signal [84]. That approach was also tested with a palladium seed and an iodine seed in an ideal environment and in beef, respectively, as a function of the seed angle relative to the incident ultrasound beam [85, 86]. It was shown from the results that the SSA algorithm outperformed B-mode images and that seeds could be detected up to a deviation of approximately 10◦ . The SSA algorithm involves many time-consuming operations, e.g. matrix computation (diagonalization, eigenvalue and eigenvector computation), eigenvalue selection, and signal reconstruction, which may cause heavy computational load.  1.2  Ultrasound Elastography  The elasticity of soft tissues depends to a large extent on their molecular building blocks and the structural organization of these blocks. Some tissues or organs are stiffer than others. Furthermore, changes in tissue stiffness are also correlated with pathologic changes. Unfortunately, current imaging modalities, such as the CT, MRI and B-mode ultrasound cannot directly measure the mechanical properties of tissues. Interest in imaging the elasticity of tissues has significantly increased in the past two decades. Elastography is a technique for quantitative imaging of strain or elastic modulus distributions in biologic tissues [98]. With this technique, external compression is typically imparted on the tissue surface with an ultrasound probe along the beam direction by a small amount (typically strain levels less than 2%). Under varying external compression on tissue, digitized ultrasound radio-frequency (RF) signals, which are the raw echoes before the formation of B-mode images in ultrasonic signal processing, are acquired and recorded, from a region of interest. An RF data frame consists of a number of scan lines, and each RF scan line is split into a number of short segments with equal length and overlap, at different depths. Tissue displacements between pre- and postcompression ultrasonic RF image pairs are estimated by various methods. The spatial derivative of the displacement field is equal to the strain in tissue. The resulting strain images or elastograms provide information about the mechanical properties of tissues. Fig. 1.5 presents the schematic diagram of elastography or strain imaging. If additional information, such as the distribution of applied force and boundary conditions, can be obtained, the distribution of elastic modulus of tissue can be reconstructed under certain assumptions using inverse-problem techniques. In medical applications, both ultrasound and MRI [71, 120] can be utilized for strain imaging. Because TRUS is intra-operatively used during the prostate brachytherapy procedure, only ultrasound elastography is studied in our research.  8  Figure 1.5: Schematic diagram of ultrasound elastography or strain imaging.  1.2.1  Tissue Motion Estimation  Tissue motion tracking plays a key role in elastography. Tissue motion can be estimated by various methods on the corresponding RF data set. In the time domain, tissue displacement can be estimated by the time domain cross-correlation (TDE) [81, 98], sum of absolute differences (SAD) or sum of squared differences (SSD) [74, 144], speckle phase tracking [97], iterative phase root seeking (PRS) [105, 106], zero-crossing tracking [124], etc. In the frequency domain, tissue displacement can be estimated by the variation in the mean-scatterer spacing [128], spectral centroid shift [69], spectral cross-correlation [139], etc. In general, the time-domain methods have the advantage of being highly precise, although analysis has shown that they are not very robust in the presence of the decorrelation between pre-compression and post-compression RF signals [141]. Phase root seeking (PRS) is a fast alternative to the conventional time domain cross-correlation approach. In contrast to the search of the maximum of the cross-correlation function, PRS searches the root (zero) of the phase of the correlation function of the corresponding analytic (complex) signals by Newton iteration [105, 106]. To estimate subsample shifts with high accuracy, interpolation is needed. The raw RF data are first converted to the analytic signals by the Hilbert transform, and demodulated to the baseband signals. Then the baseband signals are interpolated by linear interpolation. The phase root of the cross-correlation function is iteratively searched. After a small number of iterations (usually 2 or 3), the algorithm will converge and reach a certain accuracy [106]. To improve the accuracy of tissue motion estimation, logarithmic compression of the RF signals [20] was used to reduce amplitude variations caused by speckles and variations in backscatter, which cause the estimated time shift to deviate from the time shift in the center of the window. A theoretical framework, named the strain filter [141], has been developed to predict the upper bound of strain estimation in relation to the ultrasonic, mechanical and signal processing parameters. 9  Temporal stretching of the post-compression signals and multicompression averaging [8, 140] were used to reduce the decorrelation noise in elastograms from mechanical compression of tissues. Chaturvedi et al. [21] proposed a 2D local companding method, and Lorenz et al. [79] employed modified optical flow algorithm, to reduce the lateral decorrelation. Moreover, the lateral and elevational (out of plane) tissue strains resulted from the axial compression were estimated by Konofagou et al. [66] and Kallel et al. [59]. Additionally, in order to reduce the noise amplification of the gradient operation, tissue strain can be also computed by the least squares estimation (LSE) of the derivative of the displacement field [58]. A new tissue motion estimation method called time domain cross-correlation with prior estimates (TDPE) was developed by Zahiri-Azar et al. [150, 151]. Similar to the TDE method, the displacement of a tissue element (or segment) from sample time kT to sample time (k + 1)T is obtained by maximizing its cross-correlation to the RF scan line sampled at time (k + 1)T , while the prior displacement estimates of neighboring tissue elements are used as prediction to speed up this process. Because of the physical proximity, adjacent tissue elements have similar displacements when the RF image frame rate is high and strain level is relatively low. Therefore, the searching range for cross-correlation maximization can be bracketed to a very small predicted interval [151]. As this approach incorporates the prior tissue motion estimates, it works much faster and is more accurate than the conventional TDE method. To improve the accuracy, sub-sample resolution is achieved by using quadratic or cosine interpolation. The normalized cross-correlation is also monitored in real time. When the correlation coefficient drops below a certain threshold (e.g., 0.9), a “recovery search” that samples a larger search range, not just the small predicted range, is performed to obtain more reliable displacement estimates.  1.2.2  Vibro-Elastography  All elastography techniques combine mechanical excitation with measurement of the resulting tissue motion. Most strain imaging systems use quasi-static (below 10 Hz) compression [97, 98, 105, 121] and/or manual excitation. The tissue deformation is related to local stiffness of the region. The resulting strain images often have large variance, because the external compression is timevarying and operator-dependent and because the tissue motion estimation does not incorporate the information about the excitation. This kind of elastography method can only measure static behavior of the tissue but not the dynamics. In another approach, which is also called sonoelasticity imaging, an external vibration source creates high-frequency shear waves that propagate into the tissue, and the elastic properties are extracted by analyzing the wave propagation and decay [40, 75, 76, 116]. Doppler imaging was used to visualize mechanical waves to obtain tissue elasticity. Additionally, the radiation force of focused ultrasound is used to induce local tissue motion [33, 95, 143]. However, specialized high-intensity ultrasound beams are required, and the process has to be repeated for all pixel locations to obtain a full strain image. So a large image takes a long time to be acquired due to the finite tissue relaxation time. Moreover, the high intensity ultrasound beams may cause a possible hazard to patients. A dynamic ultrasound elastography method named vibro-elastography (VE) was proposed in [135, 136]. With this technique, the tissue dynamics along the axis of motion is considered as 10  a linear dynamic system, as shown in Fig. 1.6. Tissue is excited by an actuator that generates low-frequency compression waves at multiple frequencies simultaneously (typically low-pass or band-pass filtered white noise, up to 20 Hz). A number of ultrasound RF data are obtained over a short time interval, and the displacements of an arbitrary tissue element i are estimated by the TDPE method and buffered in a time sequence xi . Assuming that the response of tissue to the excitation is linear, the transfer functions Hij (ω) of the tissue displacement sequence xj with respect to the tissue displacement sequence xi viewed as an input, and the associated coherence functions Cij (ω), can be computed by well-known techniques [18, 77]. The coherence function Cij (ω) indicates the proportion of the input energy of xi at frequency ω that appears at the output energy of xj at the same ω. Its value is between 0 and 1. Therefore, the reliability of the transfer function can be evaluated by the associated coherence function. A high coherence function indicates good confidence that the system is linear and the signal-to-noise ratio is high [136]. Elastic properties of tissue such as stiffness and viscosity can be derived from the transfer functions. The magnitude of the transfer function at zero frequency |Hij (0)| represents the amount of axial static deformation of tissue element j as a response to a unit displacement of tissue i (0))| will show the relative tissue compliance between consecutive element i. Thus |Hii+1 (0)/Hi−1  tissue elements [114]. Furthermore, in order to reduce the influence of the estimation errors of tissue displacements, it is also reasonable to take an average of the magnitude of the transfer function over the range of low frequencies where it is approximately constant, instead of using the zero-frequency response |Hij (0)|. Stiffness estimation based on transfer functions relies on elastic behavior that can be averaged over an entire low frequency range and therefore it is superior to static strain images [136].  Figure 1.6: The linear dynamic tissue model in vibro-elastography, denoted by Hij , and xi and xj are the time sequences of the displacements of tissue elements i and j at different spatial locations, respectively.  11  1.2.3  Clinical Applications of Elastography  A lot of research work has been carried out in the area of clinical applications of ultrasound elastography. Garra et al. [41], Hiltawsky et al. [49] and Bercoff et al. [19] used it for breast cancer detection. The feasibility of strain imaging in intravascular applications was studied by Korte et al [26]. Cardiac elastography is another diagnostic application for detection of heart function abnormalities [17, 67, 68], such as ischemia and infarct, based on the analysis of the local cardiac muscle motion. As the prostate is a firm organ with respect to its surrounding tissues, ultrasound elastography would be a potential imaging modality for the guidance of prostate brachytherapy. In the past decade, elastographic characterizations of the canine prostate, which is widely used as a model for the human prostate, were given by Kallel et al. [60], and elastic properties of the prostate were measured before by Krouskop et al ex-vivo [70]. Souchon et al. [122, 123] used elastography to visualize the HIFU-induced lesions in vivo during and immediately after the HIFU therapy of prostate cancer. Lorenz and Pesavento et al. [80, 104] developed the first real-time elastography system in the world for prostate tumor examination. They found the specificity for cancer detection to be approximately 84% and a sensitivity of approximately 76% after the evaluation of 170 patients. Recently elastography has been used for prostate cancer detection [7, 62] and imageguided biopsy [65]. Moreover, the vibro-elastography technique has been tested for the elasticity modeling and imaging of the prostate region of a few patients [114] at the British Columbia (BC) Cancer Agency, Vancouver Centre, and the results were encouraging.  1.3  Prostate Boundary Segmentation  During and after prostate brachytherapy, the prostate may deform and swell due to the TRUS probe re-positioning, needle insertion and retraction, and edema. Thus prostate boundary segmentation also plays a significant role for prostate brachytherapy. The brachytherapy procedure depends on the accurate and real-time imaging to distinguish the prostate from the adjacent tissues. However, prostate boundary segmentation in B-mode TRUS images, whether done manually or automatically, still remains a challenging task, because of the poor contrast between the prostate gland and surrounding tissues, speckle noise, shadowing and refraction artifacts [117, 118]. Fig. 1.7 shows an example of the prostate B-mode TRUS image in the transversal scan plane. For these reasons, a fully-automatic prostate boundary segmentation method in conventional ultrasound images is very difficult to implement, and a priori knowledge of the anatomical structure of the prostate and intervention of observers are usually needed. Various prostate segmentation methods in B-mode ultrasound images have been proposed in the literature [96, 117]. Edge-based detection is most commonly used for image segmentation. The common feature of this class of methods is the use of edge detectors to identify all the edges in the image, followed by edge selection and linking to outline the boundary. Such methods usually lead to spurious boundaries in highly textured areas, and missed boundaries where the prostate is not well delineated due to the poor quality of ultrasound images [117]. Therefore, a priori knowledge of the prostate structure becomes necessary for the edge selection and linking. Chen et al. [22] 12  Figure 1.7: Example of the prostate B-mode TRUS image in the transversal scan plane, where the TRUS probe is on the bottom and is normal to the image. There are three implanted seeds close to the anterior boundary of the prostate in the image.  proposed a derivative edge detection approach and Aarnink et al. [3] used nonlinear filtering (minimum/maximum filter). Pathak et al. [103] proposed an edge-guided approach using weak membrane fitting, and Awad et al. [12] developed another method based on the Canny detector and morphological opening. In the latter two methods, an algorithm called sticks was used to enhance the boundaries and reduce speckles in the TRUS images at the same time [103]. It uses line segments (“sticks”) with the same length in different angular orientations as templates and selects the orientation at each point that is most likely to represent a line in the image. In the output image, the contrast at the edges is increased, while the noise is decreased. Most of the developed methods use deformable contour models (or active contour models, snakes). During the past several years, deformable contour models have been extensively investigated as an appealing tool for image segmentation. Pathak et al. [102], Knoll, et al. [63, 64], Ghanei et al. [42], Ladak, et al. [72], and Hodge et al. [50] presented algorithms based on deformable contour models to extract the prostate boundary in 2D or 3D B-mode ultrasound images. Nevertheless, these methods are computationally expensive and susceptible to local minima, because there is an optimization stage involved. Therefore, they are not suitable for real-time clinical applications. In addition, human intervention is often required to get a close guess to the actual boundary for the initialization of the deformable contour models. Statistical shape models have also been used for prostate contour detection. These models code the variations of a group of selected parameters that are used to describe the detected object in an observed population or training sets [117]. The parameters are estimated from the available segmented images and are used for the new images. Moreover, the results from new segmentations can be incorporated into the model, thus updating a priori information available to the system. Prater et al. [110] described a method with feed-forward neural networks. Lorenz et al. [78] utilized the Markov random field to extract the prostate boundary. Shen et al. [118] proposed a statistical shape model with Gabor filters at different orientations and scales for automatic 2D prostate boundaries detection. The authors also extended their work to segment 3D ultrasound 13  prostate images [153]. Gong et al. [45] modeled the prostate shape by fitting super-ellipses to manually segmented prostate contours. Similar to the deformable contour models, the common drawback of statistical models is the intensive computational load. Abolmaesumi and Sirouspour [5, 6] developed a new method for prostate boundary segmentation, using the interacting multiple model (IMM) estimator along with a probabilistic data association filter (PDAF). The IMM/PDAF algorithm extends the Star algorithm and spatial Kalman filter from [4], and accommodates multiple trajectory models which can differ in their noise levels and/or their structures [5]. The system can switch between the models according to a set of transition probabilities. This feature allows the automatic adjustment of the filter gain based on the boundary curvature. This approach has been tested on a series of prostate ultrasound images and the results are promising. Compared with the deformable contour models or statistical models, the IMM/PDAF method does not employ any numerical optimization technique and does not need a lot of human intervention. Thus it could be implemented in real time during the prostate brachytherapy procedure. There is a seed point inside the prostate contour that can see all the boundary points to be set manually or automatically. Although a small change in the seed point location does not significantly affect the segmentation outcome, the sensitivity of the extracted contour related to the seed point position depends on the quality of the edge points and the shape of the boundary. Based on the IMM/PDAF algorithm, Salcudean and Badiei et al. [13, 14] proposed a semiautomatic segmentation method. The main purpose is to obtain smooth, continuous contours with no sharp edges, no hourglass shapes and approximate symmetry about the median axis of the prostate, which is the requirement by the British Columbia Cancer Agency, Vancouver Centre. The approach utilizes an image warping algorithm to change the prostate contour shape from concave to convex, the IMM/PDAF algorithm to segment the warped prostate contours, ellipse fitting to find the best elliptical fit to the warped prostate, and a reverse-warping algorithm to reverse the elliptically-fit prostate boundary. The proposed method was validated using manual segmentation by an expert observer on 17 midgland, pre-operative, B-mode TRUS images. The merits of the algorithm lie in its straightforward, intuitive process but most importantly in its computational simplicity and ability to create good results despite non-ideal TRUS images [13, 14]. The algorithm was almost two times faster than manual segmentation and thus it has potential for real-time applications. Recently, Mahdavi et al. [83] extended this work from 2D to 3D prostate boundary segmentation using image warping and ellipsoid fitting. Prostate boundary segmentation in ultrasound images is still a open problem. The segmentation result in B-mode TRUS images, whether done manually or automatically, is not always reliable. Based on our study in Chapter 5, the contrast between the prostate and surrounding tissue appears higher in VE images than in B-mode images. Thus, prostate segmentation will be improved with VE techniques.  14  1.4  Research Goals and Objectives  The research background and literature review presented above have shown that the intra-operative dose feedback for prostate brachytherapy is still an open and difficult problem, in which the detection of brachytherapy seeds and the segmentation of prostate boundaries are two critical factors. The overall objective of our research work is to develop an ultrasound-based solution for accurate, reliable, and user-independent intra-operative dosimetry for prostate brachytherapy. Specific research goals for the proposed dosimetry solution are outlined as follows: • To develop, implement and validate ultrasound signal and image processing algorithms that enhance imaging and detection of brachytherapy seeds. In addition, the signal and image processing algorithm(s) should be able to perform in real time. • To design, develop and implement an ultrasound-based method that can detect implanted brachytherapy seeds in 3D ultrasound images, and to validate the approach with a prostate phantom according to a realistic treatment plan. • To design and develop algorithms and application programs for the implementation of realtime vibro-elastography and seed detection methods. Requirements for the application program include the graphical user interface, online ultrasound data access, computation and display of the desired image types, real-time data storage, and offline data analysis. • To integrate all the components into a functional prototype system, which can be validated through laboratory and/or clinical tests. • To implement and validate the tomosynthesis reconstruction algorithm for seed detection in fluoroscopic images. The tomosynthesis result can be used to validate our ultrasound-based seed detection results. In addition, related work is to be carried out with my supervisors and other graduate students in the laboratory, on the modification and improvement of the RF signal processing approaches and programs that can quickly and reliably estimate tissue displacements, on the testing and validation of the mechanical components for the proposed vibro-elastography system, and on the design and implementation of making phantoms for different research purposes.  1.5  Thesis Outline  This chapter introduced the reader to the research topic and related clinical background, and outlines the thesis goals. An overview of current literature of brachytherapy seed detection, ultrasound elastography techniques, and prostate boundary segmentation in ultrasound images, have been presented in this Chapter. In Chapter 2, the candidate methods that we have tried for enhancing the visibility and imaging of brachytherapy seeds in 2D ultrasound are described with phantom experiment results and analysis. 15  Chapter 3 proposes a seed detection solution using 3D ultrasound and using a priori knowledge of the planned needle insertion entries and seed placements, and presents a validation result using a prostate phantom with realistic implant according to a plan generated by a medical physicist and implanted by a radiation oncologist. The fast computation algorithms of the transfer functions and coherence functions and application programs for the implementation of real-time vibro-elastography are described in Chapter 4. A simplified version of vibro-elastography, named power strain imaging, is also proposed in this Chapter. Chapter 5 presents the design of the vibro-elastography imaging system for clinical uses, and the preliminary results of using the vibro-elastography system to image the prostate region of a number of patients at the British Columbia Cancer Agency, Vancouver Centre. Finally, a discussion of the proposed system and its possible impact, an overview of the main contributions of this dissertation and conclusions are presented in Chapter 6. Additional materials include photocopies of the ethics certificate and of the subject information and consent forms certificated and approved by the UBC Research Ethics Board, copies of the clinical data acquisition protocols used during the volume study and in the operating room pre-operatively, communication protocol between the real-time data acquisition program and the actuation system, and recipes of making agar and gelatin phantoms, in Appendices A - F, respectively.  16  Chapter 2  Enhancement of Seed Visibility in 2D Ultrasound This chapter describes the candidate ultrasound signal and image processing methods that we have evaluated to enhance the visibility and imaging of brachytherapy seeds in 2D ultrasound, compares the preliminary results of phantom experiments, and discusses why we chose the reflected power imaging method as opposed to other methods, to image the brachytherapy seeds.  2.1 2.1.1  Candidate Methods Reflected Power of RF Signals  In ultrasound imaging, high-frequency sound waves are transmitted into the medium and the reflected echoes are recorded [23]. The reflected energy is quite different from various scattering particles in the medium, because they have different acoustic impedances. In particular, metal brachytherapy seeds are strong reflectors, and there is much more wave energy reflected from seeds than from most of other scatterers in tissue. This property can be used to develop imaging methods in which seeds provide a strong contrast. One method of measuring the reflected energy is to directly monitor the average reflected power of the ultrasound RF signals. Every frame of RF data consists of a number of scan lines. With the reflected power method, each RF scan line is divided into a series of short segments with equal length and overlap between neighboring segments, at different depths along the RF scan line. Fig. 2.1 shows an example of dividing an RF scan line into N short segments. This is similar to the elastography techniques, while external compression is not required. In general, the segment length can be set equal to the radius of brachytherapy seeds. Each RF segment corresponds to a pixel in the computed reflected power image and the location of the pixel is at the center of the RF segment. The reflected power of RF signals can be estimated in the time and frequency domains.  Frequency Domain Fig. 2.2 shows the flow diagram for computing the reflected power of RF signals in the frequency domain. The power spectrum is computed by the squared modulus of the fast Fourier transform (FFT) magnitudes. Prior to the FFT, an RF data segment is multiplied by a window function, e.g. a Hamming or Hanning window, to reduce spectral leakage [99]. Zero padding is needed if the window length is less than the number of FFT points. The reflected power of an RF segment at 17  Figure 2.1: Example of dividing an RF scan line into a series of N short segments with equal length and overlap.  depth d, denoted by PA (d), is given by ∑q  k=p Pd (fk ) (N + 1)2  PA (d) =  (2.1)  where Pd (fk ) is the magnitude of the kth frequency component in the 2N -point power spectrum of the segment at depth d. Typically only half of the power spectrum is used because it is axially symmetric about the Nyquist frequency (half of the sampling rate) [99]. p and q are the indices of the start and end frequency components. One can choose the desired frequency range to compute the reflected power by adjusting the values of p and q. When p = 1 and q = N + 1 are used, the overall reflected power from 0 Hz to the Nyquist frequency is determined. When the frequency range is chosen to be equal to the primary bandwidth of the transducer, the in-band reflected power is computed. When the frequency range is above the bandwidth of transmission waves, the reflected power of the harmonic waves is obtained.  Time Domain Alternatively, the average power of an RF segment can be determined in the time domain. According to Parseval’s relation [99], the total energy in the FFT-based power spectrum is equal to the total energy in the corresponding n-point time sequence. Fig. 2.3 shows the flow diagram. In the time domain, the reflected power of an RF segment at the depth d is computed by ∑n PA (d) =  2 k=1 w(k)x(k)  n  (2.2)  where x(k) (k = 1, ..., n) are the samples in the RF segment, and w(k) are weights typically determined by a window function, e.g. a Hamming or Hanning window. The computation of the reflected power in the time domain is typically much faster than that in the frequency domain, as the FFT is not required. Thus it is suitable for real-time applications, although it loses the ability of choosing the desired range of frequency components. 18  Figure 2.2: Flow diagram of computing the reflected power in the frequency domain.  Autoregressive Modeling Parametric spectral analysis is a third candidate. Autoregressive (AR) modeling has become more and more popular, mainly because the superior frequency resolution and the reduced variance of the estimation [15]. In general, the AR model of order p of the sampled ultrasound RF signal is given by: y(n) = −  p ∑  ak y(n − k) + u(n)  (2.3)  k=1  where y(n) represents the measured data, ak ’s are the unknown model parameters to be estimated, and the input noise u(n) is assumed to be a zero-mean white noise process with variance σu2 . The model parameters can be estimated by the Burg algorithm, batch least-squares, recursive leastsquares and other methods as discussed in [15, 36, 43]. When the p + 1 unknowns (ak ’s and σu2 ) have been estimated, the power spectrum of the AR model can be evaluated as follows: PAR (f ) =  |1 +  T σu2 −j2πf kT |2 k=1 ak e  ∑p  (2.4)  where f is the normalized frequency, and T is the signal sampling period. While using high order models may provide more accurate approximations of the signal, this leads to spurious peaks in the power spectrum estimate and requires more computational time [36]. Moreover, the AR model matches the signal spectrum very closely near spectrum peaks, even if the order is low [43]. Thus, a second-order or third-order AR model can efficiently estimate the power spectrum of a signal. 19  Figure 2.3: Flow diagram of computing the reflected power in the time domain.  2.1.2  Squared Envelope of RF Signals  A second method for measuring the reflected energy is to monitor the squared envelope of the RF signal. In theory, the square of the envelope signal is proportional to the reflected energy of ultrasound waves from the reflective objects in the medium. If an interface is provided by the manufacturers, the envelope signal could be directly acquired from ultrasound machines. Otherwise, it can be derived from the RF signal by the conventional approach. First, the digitized RF signal x(k) is converted to the analytic signal x ˜(k) by x ˜(k) = x(k) + j x ˆ(k)  (2.5)  where x ˆ(k) is the Hilbert transform of the RF sequence x(k). An N -point (N is an odd number) finite impulse response (FIR) filter is typically used for the real-time implementation of the Hilbert ˜(k) can be easily obtained by transform [99]. Then the envelope xen (k) of the analytic signal x xen (k) =  √  x(k)2 + x ˆ(k)2  (2.6)  As the square of the envelope amplitude is utilized in this method, the operation of the square root is not required in practice and thus the computational cost is reduced. If the true envelope of the RF signal is accessible from ultrasound machines, one only needs to square the acquired envelope values. To be consistent with the reflected power imaging approach, 1D mean filtering (replacing each pixel value with the mean intensity of its neighborhood) is applied along each squared baseband envelope to smooth the curve. Finally, the filtered squared baseband envelopes for all the RF lines are grouped together into a squared envelope image and the image could be resampled onto a desired rectangular grid if needed.  20  2.1.3  Decompressed B-mode Images  A third candidate approach decompresses conventional B-mode images obtained from ultrasound machines. B-mode images are derived from the envelope of the RF data, after the scan conversion (resampling onto a rectangular grid), logarithmic compression, thresholding, quantification and other nonlinear operations built into ultrasound machines. However, the technical details of those operations are typically proprietary to the manufacturers and unknown to users. The purpose of decompressing B-mode images is not to fully recover the true envelope signal, but to increase the contrast between the implanted seeds and surroundings. To decompress B-mode images, one needs to invert the logarithmic compression. Other operations do not affect the dynamic range greatly. Some of them, such as resampling and thresholding, are even not invertible. For the logarithmic compression, several groups of researchers [24, 30, 61, 109] suggested a mapping of the form y = D ln x + G  (2.7)  where x is the original envelope, y is the compressed envelope, D is the dynamic range parameter and G is the linear gain parameter. By adjusting the parameters D and G, x can be rounded towards the nearest integer within a desired range of grey levels for display (e.g. 0 - 255). The inverse function of Eq. (2.7) is given by z = exp(  y−G ) D  (2.8)  In practice, B-mode images are decompressed by z = K exp(y/D)  (2.9)  where K is a scaling factor that is initialized by 1. Based on the research of Kaplan and Ma [61] D can be initially estimated by ˆ = D  √ 6V /π  (2.10)  where V is the variance of the original B-mode image. K and D can be further manually adjusted to obtain a satisfactory decompressed B-mode image with enhanced contrast between the seeds and background.  2.2 2.2.1  Image Post-processing Sector Reformatting  A TRUS probe provides two scan planes: the transverse scan plane (a sector plane normal to the long axis of the probe) and the sagittal or axial scan plane (a linear plane containing the long 21  axis of the probe). This step is required only for the transversal scan. In the transversal scan plane, all the RF scan lines are distributed in a sector view with different angles because of the convex shape of the crystal array. Therefore, a destination-oriented polar to Cartesian coordinate transform [131] is used for the sector reformatting of the reflected power lines. For each pixel (x, y) of the destination image, the corresponding polar coordinates (ρ, θ) of the source image are computed by  ρ =  √  x2 + y 2  θ = arctan(x/y)  (2.11)  where θ is defined as the angle of each scan line with respect to the vertical direction. The associated image value is derived from the values of its nearest neighbors in the appropriate source image by bilinear interpolation.  2.2.2  Image Enhancement  For all the candidate methods above-mentioned, after a result image is computed, median filtering is applied for image enhancement because a median filter can remove speckle and preserve sharp edges at the same time [46]. A small kernel size is chosen due to the small cross section of a brachytherapy seed. Then thresholding is applied to remove artifacts with low pixel intensities unrelated to the highly reflecting seeds and enhance the visualization of the seeds. The threshold values are determined by  T1 = αmean(I) T2 = βmean(I)  (2.12)  where I is a result image, T1 is a lower threshold limit to remove the background clutter, T2 is an upper limit to reduce the intensity difference between the implanted seeds, and α and β are two adjustable scaling factors (typically, α = 1 and β > α).  2.2.3  Beam Steering  Beam steering [23] changes the transmission angle of ultrasound beams with respect to the normal of the transducer. It is one of the basic functionalities provided for users with currently-used commercial ultrasound machines After needle insertion and retraction, brachytherapy seeds (especially the loose seeds) often become oriented in a nonparallel manner with respect to one another and to the TRUS probe. Seeds with their long axis perpendicular to the incident beam are more apt to be detected by TRUS. Thus beam steering is used in the sagittal scan plane to enhance the visibility of the implanted seeds whose long axes are not parallel to the long axis of the probe. A group of beam angles are used for this purpose, and they are separated by equal angle increments (e.g. 5◦ ) 22  within a desired range (e.g. ±20◦ ). Such beam angle pattern is axially symmetric about the normal to the transducer which is defined as the 0◦ beam angle. At each beam angle, except for 0◦ , after computing the reflected power, squared envelope or decompressed B-mode images, a destination oriented coordinate transformation is performed to convert an angled image to a rectangular image. For each pixel (x, y) of the destination image, the corresponding angled image coordinates are calculated by  x ˜ = x + y tan(θ) y˜ = y/ cos(θ)  (2.13)  where θ is the beam angle. Then the image value is obtained by bilinear interpolation from its nearest neighbors in the source image. Finally, all the resulting images with beam steering are compounded together.  2.3  Experiments  An Ultrasonix 500 RP PC-based ultrasound machine (Ultrasonix Medical Corp. Burnaby, BC, Canada) with a biplane TRUS transducer was used in the phantom experiments. The computer has an Intel Pentium 4 CPU running at 3.2 GHz and 1 GB of RAM. The biplane TRUS probe provides two scan planes: the transverse scan plane (a sector image normal to the long axis of the probe) and the sagittal scan plane (a rectangular image containing the long axis of the probe). The TRUS probe can be mounted on a rigid stage so that it can be manually translated along or rotated about its longitudinal axis. Both 8-bit B-mode TRUS images and 16-bit digitized RF signals can be acquired from that machine. The center frequency of the transducer was set to 5 MHz, which is consistent with currently recommended operating frequencies for TRUS-guided prostate brachytherapy. The sampling rate of RF signals along each scan line was 40 MHz. For the reflected power method, each RF scan line was split into a series of segments with length of 0.8 mm (40 RF samples) and 80% overlap. The axial resolution is approximately 0.16 mm along each scan line. Each RF segment was multiplied by a 40-point Hanning window and then the average power of the segment was computed in the time domain. For the squared envelope method, as the true envelope signal was not accessible from the machine, a 15-point FIR filter was used to implement the Hilbert transform. Each squared envelope line was subsampled every 8 samples in order to obtain the same axial resolution as in the reflected power image. The spatial resolution of the decompressed B-mode image was the same as that of the original B-mode image. Median filtering (kernel size of 3 × 3 pixels) was applied to the resulting image. Then the threshold value was calculated and utilized according to Eq. (2.12). Sagittal images contained 127 RF scan lines along the 55-mm linear crystal array, and the image depth was set to 40 mm. In this case, the frame rate of RF data acquisition was 17 Hz. Transverse images contained 142 RF scan lines evenly distributed in a sector view angle of 112◦ , and the radial image depth was set to 50 mm. The RF data frame rate was 12 Hz in that configuration. In order to compare the candidate methods, a prostate phantom has been constructed as shown 23  in Fig. 2.4. The phantom is contained in a Plexiglass box (10 × 10 × 10 cm3 ), and made from 12% gelatin gel and 2% cellulose (by mass) (G9282-B gelatin, S-5504 cellulose, Sigma-Aldrich Inc. St. Louis, MO, USA). Please refer to Appendix F for the recipe. A cylindrical hole in the phantom mimics the rectum and permits insertion of the TRUS transducer. Mounted outside the front wall of the box is a Plexiglass needle guide with a pattern of holes to make needles move approximately along a straight line. In total, 24 dummy seeds (No. 6734 non-radioactive seeds, Medi-Physics, Inc. Arlington Heights, IL, USA) were inserted in the phantom in three different transverse planes by using a brachytherapy needle (Bard BrachyStar seed implant needle, Covington, GA, USA). Although the needle guide was used, there were still some small angles between the axes of the seeds and the normal of the transverse plane, due to the needle deflection and phantom deformation.  Figure 2.4: Prostate phantom made of gelatin and used for the experiments.  To measure the seed localization error, the gelatin prostate phantom was scanned by a CT (Picker PQ5000 CT scanner, Picker International Inc.) at the BC Cancer Agency, Vancouver Centre. The CT slice index and thickness is 1.5 mm. The CT image size is 512 × 512 pixels with the in-slice spatial resolution of 0.293 × 0.293 mm2 for each pixel. Although the ultrasound and CT scans were not performed at the same time, the CT scan result were used as the reference assuming that the phantom deformation between the ultrasound and CT scans was negligible. Additionally, an acoustically transparent phantom (without cellulose) was made from 2% agar gel (A6924 agar, Sigma-Aldrich Inc.) to measure the reflected energy from two seed types (regular and echo seeds), with respect to the angle between the incident beam and the perpendicular of the seed. The regular seed (OncoSeed, Amersham model 6711) has a smooth cylindrical surface. The echo seed (EchoSeed, Amersham model 6733) has a corrugated surface and incorporates the use of alternating ridges and grooves in a regular sinusoidal pattern along the length of the seed except at the two ends [25, 132]. In the experiment, a regular seed and an echo seed were embedded on the surface of the phantom, and the phantom was put in a water tank. The TRUS probe was mounted on a precision rotation stage and put in the water as well. The long axis of the TRUS probe and the seeds were approximately parallel. When the TRUS probe was rotated by the precision rotation stage, the incident beam angle to the long axis of the seed was changed. The reflected 24  power images of individual seed types were calculated in the sagittal and transverse scan planes. In the sagittal scan plane, a square region of 1 × 1 cm2 approximately centered on each seed, was also analyzed by computing the sum of intensities (SOI) for that region of interest (ROI). Furthermore, beam steering was implemented with the interface provided by the manufacturer of the ultrasound machine. The beam angle in the sagittal scan plane can be set to an arbitrary angle using that interface.  2.4  Results  Figs. 2.5 (a)-(d) show the original B-mode TRUS images and corresponding reflected power, squared envelope, and decompressed B-mode images of the implanted seeds in the phantom in three different transverse scan planes, respectively. Figs. 2.6 (a)-(d) depict the original B-mode TRUS, decompressed B-mode, reflected power and squared envelope images of the seeds in the phantom in one sagittal scan plane. The image quality and contrast of the decompressed B-mode, reflected power and squared envelope images look similar. The normalized curves (normalized to 1.0 with individual peak values) along the center scan line containing two implanted seeds, extracted from the B-mode image and its corresponding decompressed B-mode, reflected power and squared envelope images in a transverse scan plane (the middle column in Fig. 2.5), are presented in Figs. 2.7 (a)-(d), respectively. Furthermore, to quantitatively measure the performance of the proposed methods, the contrast and contrast-to-noise ratio (CNR) between the implanted seeds and the background were computed in the B-mode, decompressed B-mode, reflected power and squared envelope images in Figs. 2.5 (a)-(d). First, the centroid of each seed in an image was taken after a global thresholding was applied to the image. Then a seed region was determined by selecting the pixels that are around the seed centroid and are within the cross-section area of the real seed. The background was obtained by subtracting all the seed regions from the images. Table 2.1 presents the computed contrast and CNR results between all the seed regions and the background in the three transverse scan planes in Figs. 2.5 (a)-(d). The contrast and CNR are defined as follows:  Contrast = µs /µb 2(µs − µb )2 CN R = σs 2 + σb 2  (2.14)  where µs and µb are the averages of the pixel intensities in the seed region and background in the images, and σs and σb are the standard deviations of the pixel intensities in the seed and background regions, respectively. The CNR definition was derived from [142]. Fig. 2.8 shows the three CT images of the seeds in the gelatin prostate phantom, corresponding to the three transverse scan planes in Fig. 2.5. Errors between the seed centroids in the transversal reflected power images and the corresponding CT images were calculated. The centroid of each seed in an image was taken after a global thresholding was applied to the image. Table 2.2 presents the minimum, maximum, average and median of the absolute errors. 25  Table 2.1: Contrast and CNR between all the seed regions and the background in the B-mode, decompressed B-mode, reflected power and squared envelope images, in the three transverse scan planes in Figs. 2.5 (a)-(d). B-mode Contrast CNR  2.16 21.82  Decompressed B-mode 36.05 22.54  Reflected power 59.48 52.37  Squared envelope 59.24 51.76  Table 2.2: Statistics of the absolute errors between the seed centroids in the transverse reflected power images and the corresponding CT images. Minimum 0.01 mm  Maximum 0.80 mm  Mean 0.30 mm  Median 0.32 mm  Table 2.3 compares the running time for computing a frame of decompressed B-mode, reflected power or squared envelope image, respectively, in the transverse and sagittal scan planes under the before-mentioned configurations, relative to the data acquisition frame period, using the Ultrasonix computer. The runtime of the squared envelope image included the Hilbert transform and envelope detection. Fig. 2.9 displays the 3D visualization results of the 24 implanted seeds in the phantom from two different view angles, using Matlab functions, based on a series of sagittal reflected power images acquired when the TRUS probe was slowly rotated around its long axis by hand.  Table 2.3: Running time for computing a frame of decompressed B-mode, reflected power or squared envelope image, respectively, relative to the data acquisition frame period, using the Ultrasonix computer. Scan plane Transverse Sagittal  Data acquisition frame period 83.3 ms 58.8 ms  Decompressed B-mode 13.3 ms 12.5 ms  Reflected power 24 ms 17 ms  Squared envelope 27 ms 19 ms  Fig. 2.10 compares the reflected power images of a regular seed and an echo seed on the surface of the acoustically transparent phantom in the sagittal scan plane, at beam angles from 0◦ to 25◦ away from the perpendicular in 5◦ increments. Fig. 2.11 compares the reflected power images of the two types of seeds on the surface of that phantom in the sagittal scan plane, at beam angles from 5◦ to 10◦ away from the perpendicular in 1◦ increments. Moreover, Fig. 2.12 (a) and (b) present the curves of individual seed types as determined by the SOI within the ROI from the sagittal reflected power images, at beam angles from 0◦ to 25◦ in 5◦ increments, and at beam angles from 0◦ to 10◦ in 1◦ increments, respectively. Additionally, Fig. 2.13 shows the reflected power images of a regular seed and an echo seed on the surface of that phantom, at beam angles 26  of 10◦ and 15◦ , after performing beam steering, in the sagittal scan plane.  2.5  Discussion  From the experimental results in Figs. 2.5 and 2.6, one can see that the seeds are not clearly visualized and sometimes appear blurred although the implanted seeds could be found in the Bmode TRUS images. As mentioned before, B-mode TRUS images derived from the true envelopes are subject to resampling, logarithmic compression, thresholding, quantification and/or other proprietary mappings. These processes affect the statistics of the envelope data, and a lot of the attractive properties derived from the envelope amplitude are lost [30, 61]. Logarithmic compression sacrifices the dynamic range for high-amplitude signals in order to expand the dynamic range for low-amplitude signals. Therefore, it reduces the contrast between the strong reflectors, such as brachytherapy seeds, and the surrounding tissues. In contrast, the implanted seeds look very clear in the reflected power, squared envelope and decompressed B-mode images, even when they are in the shadows of others. The proposed three techniques produce higher contrast between the implanted seeds and the background than the conventional B-mode TRUS, as shown in Fig. 2.7 and Table 2.1. The effects of shadows and double echoes caused by the implanted seeds are reduced as well. The proposed three methods can be also utilized to enhance the visibility of other objects, such as bone surfaces [147]. Among the three methods presented, decompressing the B-mode images is the simplest approach. For the B-mode images acquired from ultrasound machines, most technical details of ultrasound signal and image processing are proprietary to the manufacturers, e.g., the logarithmic compression could use different forms, and even if it uses Eq. (2.7) the values of parameters D and G are unknown. Although the true envelope cannot be fully recovered with Eq. (2.9), we can obtain higher contrast between the implanted seeds and the surroundings by finding an appropriate value of parameter D and performing an exponential expansion. Table 2.1 indicates that decompressing B-mode images may not evidently improve the CNR between the seed and the background, although it can enhance the contrast. It is nevertheless true that sometimes this method might not identify all the seeds because of the information loss from the original B-mode image. By contrast, the other two methods are based on filtering of the digitized RF data or use of the true envelope. Because the squared envelope is proportional to the energy of reflected ultrasound waves, the mean-filtered squared envelope can be also regarded as a measure of the average power of an RF segment. Therefore, the two techniques produce resulting images with approximately equal contrast and CNR between the seed and background. Table 2.1 shows that they produce higher contrast and higher CNR than decompressing B-mode images. For the reflected power method, choosing the length of an RF segment is a tradeoff between the desired spatial resolution, sensitivity and contrast. If the segment length is too small, it is affected a lot by speckle. If the segment length is too large, it blurs the profile of the seed and sacrifices the resolution. Typically, the segment length is set between the radius and the diameter of brachytherapy seeds. Furthermore, as the reflected power and squared envelope methods involve a smoothing process onto the raw RF or envelope signals with finer sampling intervals, they could produce smoother 27  curves than the corresponding B-mode and decompressed B-mode curves. For seed localization, the axial spatial resolution along each RF scan line of the two methods is determined by the center spacing of the RF segments or by the resampling interval of the squared envelope. If the center spacing or resampling interval is small enough, they can produce finer axial spatial resolution than that of the B-mode or decompressed B-mode images. The lateral spatial resolution of the three approaches depends on the density of scan lines. Additionally, the three algorithms can be implemented efficiently because no numerical optimization is involved. Table 2.3 illustrates that the running time of each approach is much less than the data acquisition frame period (1/12 ≈ 83.3 ms for the transverse scan, and 1/17 ≈ 58.8 ms for the sagittal scan). For the reflected power method, instead of using Eq. (2.2), the average power of each RF segment can be replaced by computing the squared p-norm or squared ∞-norm of that segment in the time domain. With xi as the value of the ith sample of an RF segment x containing n RF samples, the p-norm and ∞-norm are defined as ∥x∥p = (  n ∑  |xi |p )(1/p) (1 ≤ p < ∞)  (2.15)  i=1  ∥x∥∞ = max |xi | 1≤i≤n  (2.16)  While typically p = 2 defines the quadratic norm of x, different values of p are obviously associated with the overall signal strength within the RF segment. The average reflected power defined in Eq. (2.2) can be considered as a slightly modified version of the 2-norm. With p = 1, the p-norm approximates the absolute area under the RF signal in the segment. When p grows large, the p-norm approximates the maximum RF sample magnitude within the segment. In order to compute dosimetry for prostate brachytherapy, seed detection results from the reflected power, squared envelope or decompressed B-mode images have to be combined with prostate segmentation results in the corresponding B-mode TRUS images. Image registration between the decompressed B-mode image and original B-mode image is not required as they have the same spatial location and resolution. Image registration between the reflected power or squared envelope image and the corresponding B-mode image can be implemented easily by interpolation because the B-mode images are also derived from the same RF data but have different spatial location and resolution. Figs. 2.10 and 2.11 show that the reflected power from the regular and echo seeds varies as a function of the angle between the incident beam and the plane that is orthogonal to the seed axis, which is called seed angle in the context. The results are consistent with those of Davis et al. [25] and Tornes et al. [132]. Seed visibility also depends on the angular orientation of individual seeds in the reflected power images. As the seed angle increases, the reflected energy from the regular seed decreases and the image contrast decreases correspondingly as well. The image quality and contrast are acceptable when the seed angle is not larger than 10◦ . For a regular seed, only the two ends appear clearly in the sagittal scan plane when the seed angle is larger than 10◦ . The reflected energy from an echo seed does not decay with the angle in the same manner. 28  Its corrugated surface appears to enhance the reflection as the seed angle varies, which makes it easier for seed identification. The reflected energy from the echo seed is higher than that from the regular seed when the beam angle is larger than 10◦ . Davis et al. [25] and Tornes et al. [132] showed that most of implanted loose seeds are aligned within the range of ±30◦ from the optimal orthogonal position. When stranded seeds (Rapid STRAND, Amersham) are used, misalignment reduces (typically to within ±10◦ ) because the seeds and spacers are linked together in a stiffened Vicryl suture material and are approximately coincident with the needle insertion tracks. Comparing Fig. 2.13 with Fig. 2.10, one can see that by beam steering, the incident beam direction is adjusted to be approximately perpendicular to the seed and thus seed visibility is enhanced in the sagittal scan plane when the long axis of the seed is not parallel to that of the TRUS probe. In practice, it is very difficult and not needed to steer the beam angle to be orthogonal to each implanted seed. Instead, we could use a group of beam angles axially symmetric about the normal of the transducer and separated by angle increments of 5◦ . Thus, seed angles are limited within the range of ±5◦ for individual seeds at least in one of the beam steering result images. Nevertheless, beam steering requires additional time for data acquisition. It takes 1 to 2 seconds for the Ultrasonix machine to change the beam angle and recapture image data in a sagittal scan plane. Thus it would be time-costing to perform beam steering in a large number of sagittal scan planes at different angular positions. The phantom that was used in the experiments may not be sufficient to represent the real prostate region, and the number of implanted seeds in the phantom is also less than that in a pre-operative dosimetry plan. Therefore, a more specific prostate phantom would be helpful to validate a potential TRUS-based seed detection method against a realistic treatment plan.  29  (a)  (b)  (c)  (d)  Figure 2.5: (a) Original B-mode images, (b) reflected power images, (c) squared envelope images, and (d) decompressed B-mode images of the seeds in the phantom, in three different transversal scan planes (from the left to right), where the TRUS probe is on the bottom and normal to the scan planes.  30  (a)  (b)  (c)  (d)  Figure 2.6: (a) Original B-mode images, (b) reflected power images, (c) squared envelope images, and (d) decompressed B-mode images of the seeds in the phantom, in one sagittal scan planes, where the TRUS probe is on the bottom and parallel to the scan plane.  31  (a)  (b)  (c)  (d)  Figure 2.7: Normalized curves extracted from (a) the original B-mode TRUS image, (b) the corresponding reflected power image, (c) squared envelope image, and (d) decompressed B-mode image, in a transverse plane (the middle column of Fig. 2.5).  Figure 2.8: CT scan slices of the seeds in the gelatin phantom, in the three transverse planes.  32  Figure 2.9: 3-D visualization results of the 24 implanted seeds in the phantom from two different view angles, using Matlab functions, based on a series of sagittal reflected power images acquired when the TRUS probe was slowly rotated around its long axis by hand.  Figure 2.10: Reflected power images of a regular seed (top row) and an echo seed (bottom row) on the surface of the acoustically transparent phantom, at the beam angles from 0◦ to 25◦ away from the perpendicular in 5◦ increments, in the sagittal scan plane.  Figure 2.11: Reflected power images of a regular seed (top row) and an echo seed (bottom row) on the surface of the acoustically transparent phantom, at the beam angles from 5◦ to 10◦ away from the perpendicular in 1◦ increments, in the sagittal scan plane. 33  (a)  (b)  Figure 2.12: Reflected power curves of individual seed types as determined by the sum of intensities (SOI) from the sagittal reflected power images, (a) at beam angles from 0◦ to 25◦ in 5◦ increments and (b) at beam angles from 0◦ to 10◦ in 1◦ increments away from the perpendicular.  (a)  (b)  (c)  (d)  Figure 2.13: Reflected power images of a regular seed at the beam angles of (a) 10◦ and (b) 15◦ , and reflected power images of an echo seed at the beam angles of (c) 10◦ and (d) 15◦ , after performing beam steering, in the sagittal scan plane.  34  Chapter 3  Brachytherapy Seed Detection Using 3D Ultrasound This chapter proposes an automated seed detection solution using 3D ultrasound. A 3D reflected power image is computed from ultrasound radio-frequency (RF) signals, instead of conventional Bmode images. Then implanted seeds are segmented in local 3D search spaces that are determined by a priori knowledge, e.g. needle entry points and seed placements along each needle in the pre-operative dosimetry plan. Furthermore, needle insertion tracks are also detected in local 3D search spaces by the Hough Transform to enhance seed localization results. Phantom experiments have been performed to test the effectiveness of the proposed solution, according to a realistic treatment plan for a tissue-equivalent ultrasound prostate phantom.  3.1  3D TRUS Imaging System  A 3D TRUS imaging system has been constructed. The mechanical and electrical design and implementation were completed by Salcudean et al. The author designed and developed the programs for RF data acquisition and for signal/image processing, and validated the functionalities of the 3D TRUS imaging system. The system utilizes a Sonix RP PC-based ultrasound machine, with a dual plane linear/microconvex broadband 5-9 MHz endorectal transducer (Ultrasonix Medical Corp. Burnaby, BC, Canada). Both 8-bit B-mode images and 16-bit digitized RF signals can be acquired simultaneously from the machine. The spatial sampling rate of RF signals is 20 MHz. The biplane TRUS probe provides two scan planes: the transversal scan plane (a sector image normal to the long axis of the probe) and the sagittal scan plane (a rectangular image containing the longitudinal axis of the probe). As illustrated in Fig. 3.1, the endorectal TRUS probe is mounted into the cradle of a CIVCO EX II stepper (CIVCO Medical Solutions, Kalona, IA, USA). With the stepper and cradle, the TRUS transducer can be manually translated along the longitudinal axis of the probe. A gear mounted on the cradle engages the shaft of an encoder which allows for accurate sensing of the transducer para-sagittal angle θ. In our current set-up, the original EX II encoder was replaced with a motor with a gear head and its own encoder, allowing for automatic rotation of the transducer through an angle range of (-48◦ to 50◦ ). The 3D image coordinates (x, y, z) are defined in a left-handed coordinate system in which the x-y plane is the transversal scan plane, and the z-axis is along the longitudinal axis of the probe, as shown in Fig. 3.1 (b).  35  (a)  (b)  Figure 3.1: (a) Photograph of the 3D TRUS imaging system, and (b) schematic diagram showing the coordinates used and the rotation motor. The motor and cam mounted on the transducer cradle are used to impart radial motion to the TRUS for elastography images. This feature of the system is not used in this chapter.  3.2 3.2.1  3D Reflected Power Imaging 2D Reflected Power Imaging  The 2D reflected power images are computed from ultrasound RF signals using the time-domain method described in Chapter 2.  3.2.2  3D Image Reconstruction  The 3D reflected power image is reconstructed from a series of 2D sagittal reflected power images, when the probe is rotated around the z-axis. A destination-oriented approach [35, 131] is utilized here: for each voxel (x, y, z) of the destination image, the corresponding polar coordinates (ρ, θ,  36  z) of the source image are determined by (ρ, θ, z) = (  √  x2 + y 2 , arctan(x/y), z)  (3.1)  where ρ is the distance from the voxel to the probe center, and θ is defined as the angle of a sagittal slice with respect to the vertical or the y-axis, as shown in Fig. 3.1 (b). The associated voxel value is derived from the values of its nearest neighbors between the appropriate source slices by bilinear interpolation [99]. Additionally, the resulting 3D image is smoothed by a 3D weighted mean filter [46]. As mentioned before, brachytherapy seeds are cylindrical in shape. Therefore, to reduce computational costs, the 3D kernel is set to be a cuboid. The length of the cuboid is approximately equal to that of the seed. The side length of the cross section of the cuboid is approximately equal to the diameter of the seed.  3.3  Seed Detection  Seed detection is performed locally in the 3D reflected power image. Fig. 3.2 illustrates the flow diagram of the seed detection method step by step. The details are described in the following context.  3.3.1  Volume Cropping  The goal of volume cropping is to decrease the search space for seed localization. The entry point for each needle insertion can be obtained from the pre-operative dosimetry plan, and all the seeds inserted by one needle entry will stay close to the needle insertion track. This a priori knowledge can be used to reduce the search space, from the entire 3D volume to a local volume of interest (VOI) around the needle entry. The search space can be cropped to a cylindrical volume of radius r around the needle entry point. In our experiments r = 3.5 mm is used, as the minimum distance between two needle entry points is 5.0 mm when a conventional needle guide is used. Thresholding In every local search space, candidate seed voxels are determined by thresholding as follows:  If I(x, y, z) > T Then  I(x, y, z) = 1  Else  I(x, y, z) = 0  where I(x, y, z) is a voxel inside the local search space or VOI, and T is the threshold value which is given by T = αTotsu  (3.2) 37  Figure 3.2: Flow diagram of seed detection steps.  38  where Totsu is the threshold value using Otsu’s method [101], and α is a scaling factor (initialized by 0.8 in our experiment). By thresholding, the 3D gray-scale image is converted to a 3D binary image, in which all the non-zero voxels are candidate seed voxels.  3.3.2  Seed Voxel Grouping  To localize the implanted seeds, seed candidate voxels need to be grouped into a series of seed voxel clusters. This is accomplished by first finding connected components in the local search space or VOI and then calculating the distance between the connected components in the search space. Ideally, every connected component is supposed to be a seed voxel cluster. If the distance between two connected components is less than a threshold value (this value is set to 2.0 mm in our experiments), the two connected components in the VOI will be merged as one seed voxel cluster. Furthermore, the length of each connected component is also compared with twice the length of brachytherapy seeds. Thus if the length of a connected component is greater than 10 mm, the connected component will be separated into two seed voxel clusters.  3.3.3  Seed Cluster Squeezing  Due to the sector-distribution of scan lines in the transversal scan plane, the image resolution varies with the distance from a point to the probe center. The tangential image resolution decreases with distance to the probe center, and causes wider seed profiles to be seen in the transversal (x-y) plane. To normalize the seed profiles, the size of each seed voxel cluster is squeezed tangentially to an area of approximately 1.0 × 1.0 mm2 around the seed cluster center in the transversal plane.  3.3.4  Non-seed Artifact Removal  In addition to the real seeds, some spurious seeds may also be present among the various voxel clusters. As brachytherapy seeds are small cylinders with known length (4.5 mm), the spurious seeds can be removed by checking the length of the seed voxel clusters. We assume that the length of the seeds in the 3D reflected power image ranges from 2.0 to 5.0 mm, because of seed orientations. Thus, any seed cluster with length smaller than 2.0 mm will not be regarded as a candidate seed and will be removed. Furthermore, to deal with the detected seed clusters close to the boundary of the local search volume, the volume of a seed cluster is checked as well. If the volume of a seed cluster is smaller than a volume threshold, then the seed cluster is removed. The volume threshold we used is 1.0 × 1.0 × 2.0 mm3 .  3.3.5  Needle Insertion Track Detection  In pre-operative plans for brachytherapy, needle insertion tracks are straight lines that are parallel to the probe (or the z-axis) and to one another. However, in practice, due to prostate deformation and needle bending, needle insertion tracks are not always straight lines that are parallel to the probe. However, the curvature of needle insertion tracks is typically small and the orientation from the z-axis is also limited (typically less than ±100 ). Therefore, we assume that needle tracks are approximately straight lines having small orientation angles with respect to the z-axis in 3D 39  spaces. All the seeds inserted by one needle entry stay close to the needle insertion track. Thus line detection methods can be used to detect needle or seed insertion tracks. As we show later, the needle track detection results will help with seed localization. Equation of Lines in 3D Lines in 3D spaces can be represented by the symmetric form x − x0 y − y0 z − z0 = = a b c  (3.3)  where a, b, c are components of a vector parallel to the line, and (x0 ,y0 ,z0 ) are the coordinates of any point on the line. In total there are six parameters involved in Eq. (3.3). To reduce the number of parameters representing lines in 3D, the azimuth and elevation angles of the lines and the intersection point with the z=0 plane can be used. The 3D line is expressed by the parametric equation  x = x0 + z tan(θ) √ y = y0 + z 2 + (x − x0 )2 tan(ϕ) z = z  (3.4)  where θ is the azimuth angle, ϕ is the elevation angle, and (x0 ,y0 ) are the intersection point with the z=0 plane. Fig. 3.3 shows the schematic diagram of the representation of lines in 3D space.  Figure 3.3: Schematic diagram of the representation of lines in 3D spaces.  Hough Transform The Hough transform can be used to isolate features of a particular shape within an image and has been used to detect features such as needles in ultrasound images [31]. The purpose of the 40  technique is to find imperfect instances of objects within a certain class of shapes by a voting procedure [54]. This voting procedure is carried out in a parameter space, from which object candidates are obtained as local maxima in a so-called accumulator or votes space that is explicitly constructed by the algorithm for computing the Hough transform [54]. Because it requires that the desired features be specified in some parametric forms, the classical Hough transform is most commonly used for the detection of regular curves such as lines, circles, and ellipses [46]. The Hough transform is used here to detect needle insertion tracks in the cropped 3D spaces. As Eq. (3.4) is used to represent needle insertion tracks in 3D space, the votes space of the Hough transform consists of four parameters: θ, ϕ, x0 , and y0 . For a detected seed voxel (x, y, z) in each local search volume, we calculate x0 and y0 for each candidate θ and ϕ, and then increment the votes bin at (θ, ϕ, x0 , y0 ) in the 4D votes space by 1. After the voting procedure, thresholding is applied to remove false votes peaks, and then morphological dilation [46] is performed to merge adjacent votes peaks in the 4D Hough votes space. Finally, the votes bin with the highest Hough votes is chosen to represent the needle insertion track in the local search volume.  3.3.6  Seed Matching  The goal of this step is to determine both the number of seed candidates and their location along the detected needle insertion track. Following the Hough transform, the detected needle tracks are used to enhance the seed detection result. We call this process seed matching. In addition to the needle entries, the a priori knowledge of the seed placements along each needle in the pre-operative plan is used. The planned seed placements are represented by a sequence of 1’s and 0’s, where 1’s denotes seeds. The spatial resolution of the sequence is equal to the resolution along the z-axis. In a local search volume, all the candidate seed voxels in each transversal (x-y) plane are projected onto the detected needle insertion track by v(xt , yt , zt ) =  ∑  I(xi , yi , zt )  (3.5)  i  where xt , yt , zt are the coordinates of a voxel along the detected needle track, and I(xi , yi , zt ) is the candidate seed voxel in the z = zt plane in the VOI. After the projection, in each local search volume, a normalized cross-correlation is performed between the projection and the planned seed placements along the detected needle track. The maximal correlation coefficient means the best match between the projection and the planned seed placements, and the lag index associated with it is used to represent the offset between the planned and actual seed locations.  3.3.7  Result Checking against the Pre-Plan  The goal of this step is to determine a new threshold level to be used in the local search spaces of Step 2). After the best match is obtained, the number of detected seed projections is checked against the number of the seeds in the preoperative plan. There are three possible circumstances: 41  (a) If the detected seed number is equal to the planned seed number, it is concluded that all the seeds in the current search volume have been identified and there are no false positive seeds. It will exit the iteration and continue to the last step - seed localization, for the current needle entry and local search volume. (b) If the detected seed number is greater than the planned seed number, it is concluded that there must be spurious or false positive seeds in the detection result. (c) If the detected seed number is less than the planned seed number, it is concluded that there must be seeds missing in the detection result. For circumstances (b) and (c), if the correlation coefficient is high enough (greater than 0.6 in our experiments), the iteration will stop and the process will move to the last two steps of the algorithm. Otherwise, it will go back to Step 2), where the thresholding parameter α will be adjusted. For circumstance (b), α is to increased by 0.1 and for circumstance (c) α is to be decreased by 0.1.  3.3.8  False Positive Removal and/or Missing Seed Identification  For circumstance (b), if a candidate seed does not have a match in the preoperative plan, the candidate seed is considered as a spurious or false positive seed and is removed from the detection result. For circumstance (c), if a planned seed cannot find a match in the current seed detection result, a candidate seed is added to the result approximately at the planned location along the detected needle insertion track.  3.3.9  Seed Localization  After all the candidate seeds have been detected, the center of each seed is calculated by the average coordinates of the voxels within the seed cluster: ∑n x ¯=  k=1 xk  n  ∑n , y¯ =  k=1 yk  n  ∑n , z¯ =  k=1 zk  n  (3.6)  where (¯ x, y¯, z¯) are the coordinates of the seed center.  3.4  Fluoroscopic Tomosynthesis  In addition to TRUS, X-ray fluoroscopy is commonly used in prostate brachytherapy to visualize the implanted sources intra-operatively. It produces higher contrast and higher signal-to-noise ratio between metal seeds and soft tissue than ultrasound. Fluoroscopic tomosynthesis [138] was proposed to reconstruct the 3D seed distribution from the seed segmentation results in five or more 2D fluoroscopic images that are acquired from different perspectives when the C-arm is rotated (as illustrated in Fig. 1.4). In order to validate the seed detection result from the 3D reflected power image, we developed an approach for seed segmentation in 2D fluoroscopic images, and implemented the back projection algorithm for the fluoroscopic tomosynthesis. 42  For the GE OEC 9800 C-arm (General Electric Company, NY, USA) that is used at the BC Cancer Agency, Vancouver Center, the diameter of the image detector is 9 inches, and the fluoroscopic image size is 980×980 pixels. The fluoroscopic images are stored in the Digital Imaging and Communications in Medicine (DICOM) format. The DICOM file header contains information about the rotation angle of the C-arm, however, the angular resolution is 1◦ and not sufficiently high. To obtain higher resolution (0.1◦ ), a digital protractor is mounted onto the C-arm to measure rotation angles when the C-arm is rotated.  3.4.1  Seed Segmentation in 2D Fluoroscopic Images  Seeds are segmented in 2D fluoroscopic images using a method that is derived from [37]. It consists of the following steps: (1) Selecting and opening a DICOM fluoroscopic image; (2) Manually choosing a region of interest (ROI) that is small and contains all the seeds, by clicking the top-left and bottom-right corners of the ROI, to avoid the TRUS probe and to reduce the computational cost; (3) Median filtering the cropped fluoroscopic image to remove noise and preserve sharp edges (the kernel size is smaller than the seed size); (4) Reversing the gray-scale of the filtered image to make the seeds bright and to make the background dark; (5) Applying local thresholding to the image with a small-size kernel to segment seeds and obtain a binary resulting image by T = mean(I) + α(max(I) − mean(I))  (3.7)  where I is the cropped image within the ROI, and α is a scaling factor between 0 and 1. The kernel size is set to approximately double of the seed size; (6) Applying morphological operation - image dilation - to the binary image. Image dilation is mainly used to “grow” the size of the segmented seeds so that it is more likely for the back-projected lines from the segmented seeds to intersect.  3.4.2  Back Projection  Please refer to [138] for the back projection algorithm details. In our implementation, five fluoroscopic images at approximately -10◦ , -5◦ , 0◦ , 5◦ and 10◦ relative to the vertical (0◦ ) are taken. Different from [138], the GE OEC 9800 C-arm is not iso-centric in its orbital rotation plane. Fig. 3.4 shows the schematic diagram of the GE OEC 9800 C-arm geometry. The point at which the central axis of the X-ray beam intercepts its orthogonal C-arm diameter is defined as the beam axis center (BAC). The central axis of the X-ray beam is offset from the C-arm orbital rotation center (ORC) by an amount l. The distance from the X-ray source to the image detector is defined as d. The distance from the X-ray source to the beam axis center is defined as r. As a result, Eq. 43  (1) in [138] that calculates the corresponding local coordinates (u0 , v0 ) in each projection image for a given voxel located at (x0 , y0 , z0 ) becomes:  u0 = v0 =  d(x0 sin α + y0 cos α − l) r − x0 sin α + y0 cos α dz0 r − x0 sin α + y0 cos α  (3.8)  where the definitions of (x0 , y0 , z0 ), (u0 , v0 ), α and d are the same as in [138]. For an iso-centric C-arm, the BAC is coincident with the ORC.  Figure 3.4: Schematic diagram of the geometry of the GE OEC 9800 C-arm, where the definitions of (x, y, z), (u, v) and d are the same as in [138].  A destination-oriented approach is utilized for back projection: for each voxel (x0 , y0 , z0 ) in the volume of interest to be reconstructed, we find the projected pixel (u0 , v0 ) in each projection image by Eq. (3.8). As (u0 , v0 ) are typically not integers, the image value of voxel (x0 , y0 , z0 ) is derived by nearest-neighbor interpolation. The final voxel value is the sum of the back projection results from all the fluoroscopic images. If there are five fluoroscopic images used, voxels with value of 5 will be considered as the implanted seeds. Not only the centroid but also the orientation of seeds can be reconstructed from the back projection of the 2D seed segmentation results. Once the volume of interest is reconstructed, if the number of detected seeds is less than the expected number, the back projection will be repeated after performing image dilation on the seed segmentation results with a larger dilation radius.  44  3.4.3  Phantom Validation  The tomosynthesis reconstruction approach was validated using a precisely fabricated seed phantom (please refer to [56]) provided by the John Hopkins University. There are 100 dummy seeds precisely embedded in the phantom at known positions. The seed phantom was scanned by the GE OEC 9800 C-arm at the BC Cancer Agency, Vancouver Centre. Fig. 3.5 shows the original fluoroscopic image of the seed phantom at 0◦ and the seed segmentation result in a cropped region of interest. Figs. 3.6 and 3.7 show the actual seed distribution in the seed phantom and the detected seed distribution by the tomosynthesis reconstruction, respectively. All the embedded seeds could be identified by the fluoroscopic tomosynthesis. Table 3.1 presents the statistics of the seed localization errors of the tomosynthesis reconstruction. The coordinate frame is defined the same as in Fig. 3.1 (b). The root-mean-square (RMS) errors along the x-, y- and z-axis and in the x-y plane are less than 0.5 mm. The error along the y-axis is greater than those along the the x- and z-axis.  (a)  (b)  Figure 3.5: (a) Original fluoroscopic image of the seed phantom at 0◦ , and (b) the seed segmentation result in a cropped region of interest.  Table 3.1: Statistics of the seed localization errors by the tomosynthesis reconstruction. The coordinate frame is defined as the same in Fig. 3.1 (b).  along x-axis along y-axis in x-y plane along z-axis  rms error (mm) 0.17 0.46 0.49 0.20  standard deviation of errors (mm) 0.17 0.46 0.28 0.20  max absolute error (mm) 0.38 1.10 1.16 0.60  45  Figure 3.6: Actual seed distribution in the precisely fabricated seed phantom.  Figure 3.7: Detected seed distribution in the precisely fabricated seed phantom by the fluoroscopic tomosynthesis.  46  3.5  Experiments  A CIRS Model 053 tissue-equivalent ultrasound prostate phantom (CIRS Inc., Norfolk, VA, USA) was used in the experiment. The rectangular volume enclosing the “prostate” is 5.0×4.5×4.0 cm3 . A cylindrical hole through the phantom mimics the rectum and permits insertion of the endorectal TRUS probe. A pre-operative dosimetry plan was created by a medical physicist at the BC Cancer Agency, Vancouver Centre, Canada. Seed implantation was performed by a radiation oncologist at the BC Cancer Agency, according to the pre-operative dosimetry plan. In total, 136 dummy seeds along 26 needle entries were inserted to the CIRS phantom. Fig. 3.8 shows the CIRS Model 053 tissue-equivalent ultrasound prostate phantom with the implanted seeds.  Figure 3.8: CIRS Model 053 tissue-equivalent ultrasound prostate phantom with the implanted seeds.  Fig. 3.9 shows the distribution of the needle entries in the preoperative dosimetry plan (or called the pre-plan in the context) for the CIRS phantom, in which the entry point of needle 1 was moved to c4.5 during the implant procedure. Needle indices are marked on the top-left corner of every needle entry. Following the implant, the phantom was de-gassed in water in a vacuum chamber to remove air from the needle insertion tracks. The phantom with the implanted seeds was scanned by the Sonix ultrasound machine. The sagittal scan was activated for the phantom scan and the width of the sagittal scan plane was 55 mm. Sagittal RF data frames were acquired approximately every 0.4◦ from −48◦ to 50◦ . Every RF date frame contained 128 RF scan lines. The center frequency of the transducer was set to 5 MHz, and the data capturing depth was set to 60 mm. The RF data frame rate was 20 Hz. For the reflected power imaging, each RF scan line was divided into 150 segments with length of approximately 0.8 mm (20 RF samples) and with 50% segment overlap. Thus the axial resolution of the sagittal reflected image was 0.4 mm and the lateral resolution was 55/128 ≈ 0.43 mm. To measure the seed localization error with the proposed method, the CIRS phantom was scanned by a CT (Picker PQ5000 CT scanner, Picker International Inc.) at the BC Cancer 47  Figure 3.9: Distribution of the needle entries in the preoperative dosimetry plan for the CIRS phantom. Needle indices are marked on the top-left corner of every needle entry. The number of seeds by each needle is marked inside the contour of each needle entry. The entry point of needle 1 was moved to c4.5 during the implant procedure.  Agency, Vancouver Centre. The CT slice index and thickness is 1.5 mm. The CT image size is 512 × 512 pixels with the spatial resolution of 0.293 × 0.293 mm2 for each pixel. Although the ultrasound and CT scans were not performed at the same time, the CT scan result were used as the reference assuming that the phantom deformation between the ultrasound and CT scans was negligible. In addition, the CIRS phantom was also scanned by X-ray fluoroscopy (GE OEC 9800 C-arm, GE Company) at the BC Cancer Agency, Vancouver Centre. The implanted seeds were detected by the tomosynthesis-based localization approach [138]. The fluoroscopic scan can provide higher spatial resolution (0.233 × 0.233 mm2 for each pixel) than the CT scan. In order to test the proposed seed detection method in clinical environment, three patients were scanned intra-operatively by the 3D TRUS imaging system at the BC Cancer Agency, Vancouver Centre. RF data were collected in a number of sagittal scan planes, approximately every 0.4◦ from −48◦ to 50◦ when the TRUS transducer was automatically rotated in the patients’ rectum, immediately after the implant in the operating room. The parameters used for data acquisition and image computation were the same as the phantom experiment. Fig. 3.10 shows the distribution of the needle entries in the preoperative treatment plan for one of the patients. There are 122 seeds along 22 needle entries in the pre-plan for that patient.  48  Figure 3.10: Distribution of the needle entries in the preoperative dosimetry plan for a patient. Needle indices are marked on the top-left corner of every needle entry. The number of seeds by each needle is marked inside the contour of each needle entry.  3.6  Results  Fig. 3.11 shows the seed detection result in the CIRS prostate phantom from the CT scan. All the seeds were detected successfully by the CT scan. Fig. 3.12 shows the seed detection result in the phantom from the fluoroscopic tomosynthesis. All the seeds were detected successfully by the tomosynthesis reconstruction as well. Two of the implanted seeds (the fourth and fifth seeds along needle 17) were next to each other due to seed migration by needle forces. They cannot be separated from the CT scan or fluoroscopic tomosynthesis. Figs. 3.13 (a) and (b) compare the original B-mode TRUS image and reflected power image of the seeds in the CIRS prostate phantom in one sagittal scan plane. Figs. 3.14 -3.19 show the detection results of the implanted seeds and/or needle insertion tracks in the 3D reflected power image, after the threshold, seed voxel grouping, seed cluster squeezing, non-seed artifact removal, needle insertion track detection, false positive removal and missing seed identification, respectively. All the 26 needle insertion tracks were successfully detected. All the implanted seeds were finally identified after the false positive removal and missing seed identification based on the preoperative plan. Errors between the detected seed centers in the 3D ultrasound and CT (or fluoroscopic) images were calculated for the analysis of the seed localization accuracy. Along the x-, y-, and z-axis, the errors were computed by summing up the differences between the coordinates of the seed centers. In the x-y plane, the error was measured as the Euclidean distance between the seed centers. Table 3.2 presents the statistics of the seed localization errors between the 3D ultrasound and CT scans. The statistics of the seed localization errors between the 3D ultrasound and fluoroscopic 49  Figure 3.11: Seed detection result from the 3D CT scan.  Figure 3.12: Seed detection result from the fluoroscopic tomosynthesis.  50  (a)  (b)  Figure 3.13: (a) Original B-mode TRUS image and (b) reflected power image, of the seeds in the CIRS phantom in one sagittal scan plane, where the TRUS probe is on the bottom and parallel to the scan plane.  tomosynthesis is given in Table 3.3. Table 3.2: Statistics of the seed localization errors between the 3D ultrasound and CT scans.  along x-axis along y-axis in x-y plane along z-axis  rms error (mm) 0.86 0.76 1.15 1.13  standard deviation of errors (mm) 0.87 0.76 0.45 1.14  max absolute error (mm) 1.70 1.71 2.21 3.41  Additionally, Figs. 3.20 and 3.21 show the dose curves (at the prescription dose 144.0 Gy) computed from the ultrasound, CT and fluoroscopic seed detection results, using the algorithm defined in [113], from z = 5 mm to z = 50 mm, at intervals of 5 mm. The statistics of the dose curve errors between the 3D ultrasound and CT scan, and between the 3D ultrasound and fluoroscopic tomosynthesis are presented in Table 3.4. Figs. 3.22 and 3.23 show the preliminary detection results of the implanted seeds and the needle insertion tracks, respectively, in the 3D reflected power image of the patient. Unfortunately, less than half of the implanted seeds were detected successfully by the proposed approach. The implanted seeds close to the TRUS probe are more likely to be identified. Most of the missing seeds are far from the transducer and close to the anterior or lateral boundary of the prostate. There were also false positives in the seed detection result. Needle insertion tracks 1 - 8, 12, and 16 - 22 were detected close to their pre-planned locations. Needle insertion tracks 9 and 13, 10 and 14, and 11 and 15 are very close to each other, and also have relatively large deviation angles from the pre-plan. Therefore, they could not be identified and distinguished from each other. This patient’s result is typical of all the three patients that we have scanned at the BC Cancer Agency, Vancouver Centre. 51  Figure 3.14: Detection result of the implanted seeds in the 3D reflected power image, after thresholding.  Figure 3.15: Detection result of the implanted seeds in the 3D reflected power image, after seed voxel grouping.  52  Figure 3.16: Detection result of the implanted seeds in the 3D reflected power image, after seed cluster squeezing.  Figure 3.17: Detection result of the implanted seeds in the 3D reflected power image, after non-seed artifact removal.  53  Figure 3.18: Detection result of needle insertion tracks (red lines) by the Hough transform in 3D spaces.  Figure 3.19: Detection result of the implanted seeds in the 3D reflected power image, after false positive removal and missing seed identification.  54  CT Fluoro US  300  250  250  200  200  150  150  100  100  50  50  50  100  150  200  (a)  250  300  350  400  CT Fluoro US  300  50  250  200  200  150  150  100  100  50  50  100  150  200  (c)  250  300  350  400  CT Fluoro US  300  200  200  150  150  100  100  50  50  150  200  (e)  250  300  350  400  200  (b)  250  300  350  100  150  200  (d)  250  300  350  400  400  CT Fluoro US  300  250  100  150  CT Fluoro US  50  250  50  100  300  250  50  CT Fluoro US  300  50  100  150  200  (f)  250  300  350  400  Figure 3.20: Dose curves (at the prescription dose 144.0 Gy) computed from the ultrasound, CT and fluoroscopic seed detection results, from (a) z = 5 mm to (f) z = 30 mm, at intervals of 5 mm.  55  CT Fluoro US  300  250  250  200  200  150  150  100  100  50  50  50  100  150  200  (a)  250  300  350  400  CT Fluoro US  300  50  250  200  200  150  150  100  100  50  50  100  150  200  (c)  250  300  350  400  100  150  200  (b)  250  300  350  400  CT Fluoro US  300  250  50  CT Fluoro US  300  50  100  150  200  (d)  250  300  350  400  Figure 3.21: Dose curves (at the prescription dose 144.0 Gy) computed from the ultrasound, CT and fluoroscopic seed detection results, from (a) z = 35 mm to (d) z = 50 mm, at intervals of 5 mm.  56  Table 3.3: Statistics of the seed localization errors between the 3D ultrasound and fluoroscopic tomosynthesis.  along x-axis along y-axis in x-y plane along z-axis  rms error (mm) 0.75 0.81 1.11 1.18  standard deviation of errors (mm) 0.75 0.82 0.50 1.19  max absolute error (mm) 1.81 1.86 2.64 3.60  Table 3.4: Statistics of the dose curve errors between the 3D ultrasound and CT scan, and between the 3D ultrasound and fluoroscopic tomosynthesis.  TRUS vs CT TRUS vs Fluoro  Mean absolute distance 1.29 mm 1.00 mm  Max absolute distance 4.45 mm 3.70 mm  Mean percent area error 2.61% 2.27%  Max percent area error 4.81% 3.89%  Figure 3.22: Detection result of the implanted seeds in the 3D reflected power image of the patient.  3.7  Discussion  Conventional B-mode TRUS images are subjected to scan conversion (resampling onto a rectangular grid), logarithmic compression, thresholding and other nonlinear mappings proprietary to manufacturers. Logarithmic compression sacrifices the dynamic range for high-amplitude signals in order to expand the dynamic range for low-amplitude signals. As a result, it reduces the contrast between the strong reflectors, such as brachytherapy seeds, and the surrounding tissues. By 57  Figure 3.23: Detection result of the needle insertion tracks (red lines) by the Hough transform in the 3D reflected power image of the patient.  contrast, the reflected power imaging method directly monitors the reflected energy of raw RF signals and does not logarithmically compress the RF signals. Thus it can enhance the contrast between the implanted seeds and background. In practice, an implanted seed may not be totally placed in one sagittal scan plane and 2D reflected power images cannot visualize the distribution of the implanted seeds. Therefore, a 3D reflected power image reconstructed from a series of 2D slices would provide more information for seed detection. Seed detection in a 3D reflected power image is also influenced by many factors: the distance from an individual seed to the TRUS probe, seed positions relative to one another, seed movement caused by needle retraction, other strong reflectors in the medium (such as calcification or air bubbles), etc. For the CIRS phantom, there are a few air bubbles remaining in the phantom during seed implantation, although the phantom has been degassed. The air bubbles can cause spurious or false positive seeds in the detection result. The 52 implanted seeds by needles 16, 17, and 19 - 26 (please refer to Fig. 3.9 for the needle indices) can be successfully detected even without using Step 9 - false positive removal and/or missing seed detection. These seeds are close to the probe and are not affected by shadowing from other implanted seeds. However, it is difficult to detect all the implanted seeds by needles 1 - 8, 9, 12, 15 and 18, because they are relatively far from the TRUS probe and close to the anterior or lateral boundary of the “prostate” and are interfered with by the other implanted seeds. The proposed seed detection method in the 3D reflected power image is an iterative procedure. Thresholding plays a key role in the iteration. It converts a 3D gray-scale image to a 3D binary image and produces the candidate seed voxels for the following steps. The threshold value in Eq. (3.2) is changed at each iteration by adjusting the scaling factor α. Based on the seed matching result against the preoperative dosimetry plan, α is increased by 0.1 every time if there are more seeds detected than the pre-plan, and is decreased by 0.1 every time if there are fewer seeds detected than the pre-plan. In our experiments, α is finally between 0.4 and 1.0. Nevertheless, thresholding 58  itself cannot detect all the implanted seeds correctly, especially when the needle entries are far from the TRUS probe, e.g. needle entries 1 - 8. In this case, the false positive seeds are removed from the detection result and the missing seeds are added to the result based on the a priori knowledge. Our work shows that the a priori knowledge of the pre-operative dosimetry plan can be of significant help for seed detection. First, it is used to crop each local search volume based on the planned needle entries. The local search improves both the accuracy and the efficiency of seed detection in 3D space. Secondly, it is used to determine whether all the implanted seeds are detected along each needle insertion track by checking the number of detected seeds and the seed placements and by comparing it with the pre-plan. Thirdly, it is used to remove false positive seeds and to identify missing seeds against the planned seed placements. It is very difficult, if at all possible, to detect all the implanted seeds without using the a priori knowledge. For the implanted seeds by needles 1 - 8, 9, 12, 15 and 18, there are 8 false positive seeds existing in the detection result and 3 implanted seeds missing from the detection result. They have to be removed and identified by Step 9 with the a priori knowledge of planned seed placements. The Hough transform is used for needle insertion track detection. As mentioned before, the azimuth and elevation angles of needle insertion tracks are typically small (less than ±10◦ ), and the number of intersection points on the z = 0 plane is also limited within a local volume of interest. Therefore, although the Hough votes space is 4D, the voting procedure for the Hough transform does not lead to high computational cost. The detected needle insertion tracks are used for matching the detected seeds with the planned seed placements. Additionally, the detected needle insertion tracks could also be used to speed up the registration between 3D ultrasound and fluoroscopic images. From Tables 3.2 and 3.3, one can see that for seed localization in the 3D reflected power image, the root-mean-square (RMS) errors along the x- and y-axis are less than 1.0 mm relative to the locations computed from the CT scan and from the fluoroscopic tomosynthesis. The RMS errors in the transversal (x-y) plane and along the z-axis are less than 1.5 mm. The error along the z-axis is larger than that along the x- and y-axis, primarily due to the seed orientations and reflected power variation from the seed surfaces [25]. Moreover, as presented in Table 3.4, the mean absolute distance of the dose curve errors is less than 1.5 mm and the mean percent area error of the dose curves is less than 3%, between the 3D ultrasound and CT scan and between the 3D ultrasound and fluoroscopic tomosynthesis. There are no specific guidelines for seed localization in interstitial brachytherapy. At the BC Cancer Agency, Vancouver Centre, localization errors of 1.5 mm for the transverse and 2.5 mm for the longitudinal axes are considered to be acceptable. Because in clinical practice, dose metrics such as the D90 and V100 and dose conformity are factors that enter treatment planning decisions, the location of isodose curves relative to the prostate may be a useful metric, which is why we have computed and displayed the prescription isodose errors. While successful in phantoms, the proposed seed detection method remains to be proven in patient studies. The preliminary patient results are not as good as the phantom results, and indicate that seed detection in patients is more complex and more difficult than seed detection in phantoms due to the following reasons: (1) Needle insertion tracks may have large deviation from the pre-plan although they are 59  approximately along straight lines, and even implanted seeds may not be along a straight needle insertion track, due to prostate swelling and deformation. Such large deviation or nonlinearity may cause two adjacent needle insertion tracks very close to each other and not in parallel with each other, and may even cause crossover between them (2) Seed orientation may be larger than that in phantoms due to prostate swelling and deformation, and may cause the loss of signal strength of seeds (3) Other strong reflectors in the prostate region, primarily the calcifications, may cause false positives (4) Shadowing in the prostate region may cause the loss of signal strength of seeds As mentioned in Chapter 2, beam steering can enhance the visibility of seeds with large orientation relative to the the incident ultrasound beam. However, it would be time-costing for volume data acquisition. Typically it takes 1 to 2 seconds for the Sonix machine to change the beam angle and recapture image data. In our patient studies, RF data were captured approximately every 0.4◦ from −48◦ to 50◦ . Thus, there are in total 245 sagittal scan planes. Assuming that it only takes 1 second to complete beam steering and RF data capture and that three different beam angles are used in each sagittal scan plane, the total time of data collection would be over 12 minutes, which is not practical in the operating room. Currently the ultrasound machine used during the implant in the operating room is not our Sonix RP machine, because the Sonix RP machine has not been approved for clinical use in prostate brachytherapy. Thus, we could not have collected patient data intra-operatively during the implant. A B-K ultrasound machine (B-K Medical ApS, Herlev, Denmark) is being used for prostate brachytherapy at the BC Cancer Agency, Vancouver Center. Currently the real-time access of RF signals is not provided by the B-K machine. A frame grabber with a research package will be available later for the real-time acquisition of RF data. Typically a preoperative treatment plan for prostate brachytherapy contains an even number of needle entries and these needle entries are symmetrical about the vertical, as shown in Figs 3.9 and 3.10. The implant is usually performed from the top down by a radiation oncologist at the BC Cancer Agency, Vancouver Centre, so that the seeds close to the anterior prostate boundary will not be hampered by shadowing of the seeds close to the posterior boundary and the transducer. If the frame grabber and research package for the B-K machine are ready or if the Sonix RP machine could get approved in the future, data collection can be redone every time a radiation oncologist performs 1 to 4 needle entries. In this way, only the implanted seeds by the most recent needle entries need to be detected every time, and the whole seed detection process can be implemented from the top down stage by stage. As the search volume and the number of seeds to be identified are limited every time, the probability of finding all the implanted seeds increases. In addition, a signal saturation problem with the Sonix RP machine has been found recently. The saturation causes sudden phase changes in RF signals when the received RF signal magnitude is very high. This will do harm to the reflected power computation as the acquired RF signal amplitude from the metal seeds may be lower than it true value. As a result, the contrast between the seeds and background decreases. We are expecting Ultrasonix Medical Corp. to resolve this problem. Otherwise, we may need to use other ultrasound machines that have a data acquisition 60  interface for RF signals, e.g. the B-K ultrasound machine. The signal saturation issue does not happen with the B-K machine. The image quality of the B-K machine is also higher than that of the Sonix RP machine. Overall, the experimental results show that the proposed 3D TRUS-based seed detection solution has promise in a clinical setting but is still on the way to clinical use. The preliminary patient results also indicate that needle insertion tracks that are close to the TRUS transducer but not close to one another can be detected successfully by the Hough transform. Since fluoroscopy is usually used in prostate brachytherapy as well, this will contribute to the registration between the 3D TRUS-based seed detection and fluoroscopic tomosynthesis results, by removing the requirement of fiducial markers.  61  Chapter 4  Fast Computation Algorithms and Programs for Vibro-Elastography Vibro-elastography (VE) has the primary advantages of providing more complete dynamic tissue description in terms of transfer functions (or frequency responses) and coherence functions and of producing strain images with high signal-to-noise ratio [135, 136]. The efficiency of tissue motion estimation can be greatly improved by the time-domain cross-correlation with prior estimates (TDPE) [150, 151] method. However, the estimation of transfer functions (TFs) and coherence functions (CFs) is still a significant computational load, because the Fourier transform of a sequence of tissue displacement estimates over a short time interval is involved. This chapter proposes fast transfer function estimation and coherence function estimation algorithms for the implementation of real-time vibro-elastography.  4.1  Transfer Function Estimation  The tissue model for VE is shown in Fig. 1.6, Chapter 1, where xi and xj represent the time sequences of tissue motion estimates of tissue elements i and j at different spatial locations, respectively, along an axial RF scan line. xi and xj are taken as the input and output of a linear dynamic system Hij . The transfer function Hij (ω) and associated coherence function Cij (ω) between element i to element j are computed by analyzing the power spectrum of the tissue displacements xi and xj along the axis of the compression given as  Hij (ω) = Cij (ω) =  Pij (ω) Pii (ω) |Pij (ω)|2 Pii (ω)Pjj (ω)  (4.1)  where Pii (ω) and Pjj (ω) are the power spectral densities (PSDs) of xi and xj , and Pij (ω) is the cross spectral density (CSD) between xi and xj . For tissue element i, every time a new frame of RF data is acquired, a new tissue displacement estimate is computed and then added to the time sequence xi . As time goes on, the sequence grows longer and longer, and thus the computational cost increases. It is time-consuming to compute the Fast Fourier transform (FFT) of a long time sequence directly. Furthermore, it is not feasible to save all the tissue motion estimates due to the limited memory size. By using Welch’s averaged periodogram method [99], the time sequence xi can be divided into a series of short windows with 62  equal window length and overlap. Typically the window length is a power of two for the FFT computation. The sum of PSDs and the sum of CSDs of the windows are computed and buffered. If the time interval (measured by the number of RF data frames) for TF estimation is equal to the interval between neighboring window centers, the PSD and CSD of the current data window are also available for TF estimation next time. A simple example is given in Fig. 4.1. Assuming that at time k, the PSDs and CSDs of windows A, B and C have been calculated, at time k + 1, those results can be directly reused and only the PSD and CSD of window D needs to be computed. In this way, we do not need to calculate the PSDs and CSDs of all the windows every time, as the PSDs and CSDs of the previous windows stored in the memory can be reused. This enables us to implement fast and recursive TF estimation using moving windows or using accumulative windows with a forgetting factor. Moreover, fast CF estimation can be implemented in the same way to evaluate the reliability of the transfer functions.  Figure 4.1: Example of overlapping windows used for vibro-elastography.  4.1.1  Moving Windows  In this approach, the number of windows of tissue motion estimates used for TF/CF estimation is constant. Assuming that there are N windows used, we need memory buffers to store the PSDs and CSDs of N − 1 historical windows and to store the sums of the PSDs and CSDs. Fig. 4.2 shows the flow diagram for TF/CF estimation with moving windows. Every time a new window is loaded, first the PSD and CSD of the new window are computed. Prior to the FFT, the new window is multiplied by a windowing function, e.g. a Hamming or Hanning window, to reduce spectral leakage [99]. Zero padding is needed if the window length is less than the number of FFT points. Then the PSD and CSD of the oldest window are subtracted from the PSD and CSD sums and deleted from the buffer. The PSD and CSD of the new window are saved in the buffer and added to their sums. Finally the transfer function and/or coherence function are determined with the updated sums of PSDs and CSDs. A buffer pointer is used to mark the oldest PSD and CSD in the memory to improve the computational efficiency. The PSD and CSD of the new window override those of the oldest window, and the pointer moves to the next buffer. If the pointer currently points to the last buffer, the pointer moves to the first buffer.  63  Figure 4.2: Flow diagram for the TF/CF estimation with moving windows.  4.1.2  Accumulative Windows  In another method, all the historical windows of tissue motion estimates are used for TF/CF estimation. An exponential forgetting factor λ (0 < λ < 1) is used to discount the importance of old data windows. Only the sums of PSDs and CSDs need to be stored in memory buffers. The PSDs and CSDs of individual windows are not required. Fig. 4.3 shows the flow diagram for TF/CF estimation with accumulative windows and a forgetting factor. Every time a new window is loaded, first the PSD and CSD of the new window are computed. Prior to the FFT, windowing and/or zero padding are performed as well. Then the PSD and CSD of the new window are added to the product of the forgetting factor λ and the current sums of PSDs and CSDs, respectively. As a result, the PSDs and CSDs of “old” data windows gradually lose their importance as time goes on, and the weight of the newly-computed PSD and CSD are always the highest. The smaller λ is, the faster the importance of the “old” data windows decreases. Finally the transfer function and/or coherence function are determined with the updated sums of PSDs and CSDs. As only the sums of PSDs and CSDs are stored in memory and updated every time, this approach is more efficient and saves more memory than using moving windows. 64  Figure 4.3: Flow diagram for the TF/CF estimation with accumulative windows.  4.1.3  Average of Transfer Functions  In order to reduce the influence of the estimation errors of tissue displacements, it is reasonable to take an average of the magnitude of the transfer function over the range of excitation frequencies [136]. In practice, for band-pass filtered excitation, the high cutoff frequency is limited by half of the maximum RF data frame rate of ultrasound machines, and the low cutoff frequency is close to 0 Hz. Every time the sums of PSDs and CSDs for tissue element i are updated, the average of the magnitude of the transfer function between element i to element j over the excitation frequency range is determined by ¯j = H i  ∑ j ∑ Pij (k) 1 1 | |Hi (k)| = | q−p+1 q−p+1 Pii (k) q  q  k=p  k=p  (4.2)  where p and q are the starting and ending FFT frequency components, respectively, corresponding to the excitation frequency range. After that, the least squares estimation (LSE) is applied to the result to obtain TF-based strain images (also called VE images in the context). Alternatively, an L2 norm of the differences between TFs at different spatial locations can be used to approximate tissue compliance [114].  65  4.2  Power Strain Imaging  For TF and CF estimation, both the PSDs and CSDs of the time sequences of tissue motion estimates are required. In general, the computational time for CSD estimation is a little longer than the PSD estimation time. In order to further reduce the computational load, a simplified version of vibro-elastography, named power strain imaging, is proposed here. It is based on VE techniques, in which compression is generated on tissue by an actuator excited by low-pass or band-pass filtered white noise. For tissue element i, only the PSD of the tissue displacement sequence xi needs to be computed. The PSD estimation can be implemented in the same way by using moving windows or accumulative windows. As the CSD computation is removed, the computational cost is reduced at least by half. The square root of the average power over the range of the excitation frequencies (please refer to Eq. (4.3)) is used as a measure of the tissue displacement, because the average of the power spectrum is proportional to the squared amplitude of tissue motion. Then tissue strain is determined by the LSE of the derivative of the results. The derived strain image is called the power strain image in the thesis. ∑ 1 Pii (k) q−p+1 q  P¯ii =  (4.3)  k=p  4.3  Application Programs  Two graphical user interface (GUI) based application programs have been developed for RF data analysis offline and for real-time data acquisition, respectively. The two application programs were developed with Microsoft Visual C++ 6.0 on the Windows XP platform (Microsoft Corporation, Seattle, WA).  4.3.1  RF Data Analysis Program  The RF data analysis program was designed for processing and analyzing collected RF data, displaying computed images, and saving the results images to files. The RF data are saved in binary files (*.dmp). Currently the RF data analysis program supports computation and display the following image types from the RF data: (1) Reflected power image, (2) Squared envelope image, (3) Axial displacement image, (4) Axial strain image, (5) Axial shear strain image, (6) Correlation coefficient image, (7) TF strain image, (8) CF image, and (9) Power strain image. 66  Please refer to Chapter 2 for methods (1) and (2). The C++ code for (3) axial displacement estimation using TDPE, and for estimation of (4) axial strain and (5) axial shear strain, have been provided by my colleague Mr. Zahiri-Azar [150]. Fig. 4.4 shows the main GUI of the RF data analysis program. It consists of the main frame with the menu and tool bar, the control panel on the left side, and the image client window on the right side. The control panel contains a number of buttons, checkboxes, listboxes and scrollbars for the fundamental operations, e.g. image types, color map, linear or sector view, horizontal and vertical flip, image intensity scale, temporal averaging or forgetting factor, etc.  Figure 4.4: Main GUI of the RF data analysis program, which consists of the main frame with the menu and tool bar, the control panel on the left side, and the image client window on the right side. The image properties dialog (or the advanced dialog), as shown in Fig. 4.5, can be activated by clicking the “Advanced” button on the control panel. It is used to change the parameter values to determine the dimensions of result images, to compute different types of result images, and to set the name and data type of result files. The axial resolution of the resulting images is determined by the block spacing and block length in the RF image properties dialog. Their unit is RF signal 67  samples. When the axial RF sampling rate is set to 20 MHz, 26 RF samples are approximately equivalent to 1.0 mm. Therefore, the depth of result images (in mm) is estimated by Depth =  BlockN umber × BlockSpacing 26  (4.4)  Typically, the result image depth is less than or equal to the RF data capture depth.  Figure 4.5: Advanced dialog of the RF image properties with the RF data analysis program.  The image client window displays resulting images computed with different imaging methods, and also shows related RF data information and result image properties. Fig. 4.4 shows an example of axial strain images. The display region (the rectangle in the middle) is specially tailored to fit to the dimensions of the resulting images, by using the related Ulterius software development kits (SDKs) provided by the Ultrasonix Medical Corp. Burnaby, BC, Canada. The computed result image data are also saved in binary files (*.dmp) with either float or byte data types. The result images shown in the display region can be saved in 24-bit color bitmaps. In general, the use of the RF data analysis program involves the following steps: (1) Opening the file dialog by clicking the menu or tool bar and selecting the RF data file (*.dmp) that the user wants to use. (2) Choosing the image type for computation and display from the Display Options list box on the control panel. (3) Changing the options or settings on the control panel for the color map, linear or sector view, horizontal and vertical flips, amplification, temporal averaging or forgetting factor, etc. (4) Clicking the “Advanced” button to open the RF image properties dialog to change the 68  result image dimensions, to set computational parameters for computing different types of images, and to set parameters for saving result data or images. (5) Checking the corresponding options on the control panel, if one wants to save the result data in binary files (*.dmp) or to save the result images in bitmaps. (6) Clicking the “Play” button on the tool bar to start computing and displaying result images. One can pause, resume, stop, and restart playing whenever he or she wants.  4.3.2  Real-time Data Acquisition Program  The real-time application program was designed primarily for B-mode and RF data acquisition (DAQ). Its main framework was derived and modified from the original real-time strain imaging program courtesy of Patel and Zahiri-Azar. It is installed and running on the ultrasound machines. The B-mode and RF data are retrieved from the ultrasound machines by using the related Ulterius SDKs provided by the Ultrasonix Medical Corp. Fig. 4.6 shows the main GUI of the real-time DAQ program. It also consists of the main frame with the menu and tool bar, the control panel on the left side, and the image client window on the right hand. Different types of result images (the same as listed in Section 4.3.1) can be calculated in real time and shown over the corresponding B-mode images simultaneously in the image client window. If high frame rates are preferred for saving RF data, computation and display of the result images can be limited to a small rectangular region or can be disabled (only B-mode images are shown to monitor the scan planes in real time). The RF or B-mode data are stored in memory during data acquisition and are output to binary files (*.dmp) after the data collection is stopped. Thus, a large memory size is required. This is not an issue at present. The real-time DAQ program requires 1.8 GB memory space: 1.6 GB for RF data storage and 0.2 GB for B-mode data storage, as the B-mode data size is much less than the RF data size. Moreover, the real-time DAQ program is able to communicate with the control program of the latest VE actuation system, which is described in Chapter 5, using two shared text files for the rotation angles and vibration positions of the transducer, respectively. The DAQ program keeps reading the shared text files every time a new frame of RF data is acquired. The rotation angle and vibration position values are stored in the header of every RF or B-mode data frame. These values are also shown in the image client window in real time. Additionally, the start and stop of the RF data acquisition procedure can be automated by the reading of rotation angles, when the transducer is automatically rotated by the actuation system. Please refer to Appendix E for the communication protocol between the real-time DAQ program and the control program of the VE actuation system. Typically, the use of the real-time DAQ program involves the following steps: (1) Connecting to the ultrasound server installed with the Sonix RP ultrasound machine. (2) Choosing the image type for display from the Display Options list box on the control panel. (3) Changing the options or settings on the control panel for the color map, linear or sector view, horizontal and vertical flips, amplification, temporal averaging or forgetting factor, etc. 69  Figure 4.6: Main GUI of the real-time data acquisition (DAQ) program, which consists of the main frame with the menu and tool bar, the control panel on the left side, and the image client window on the right side. (4) Changing the position and size of the VE imaging region by clicking the mouse and then dragging the rectangle, or selecting to show B-mode images only (to obtain high frame rates) on the control panel. (5) Selecting the data acquisition mode from the DAQ Mode listbox on the bottom of the control panel. Typically, “VE + Auto-rotation” is used for the VE imaging and “NonVE + Autorotation” is used for the reflected power imaging. (6) Checking/unchecking the “Save Data” checkbox to start/stop data acquisition when the probe is manually rotated (This step is not needed when the transducer is rotated automatically by the VE actuation system).  70  4.4  Experiments  An Ultrasonix 500 RP PC-based ultrasound machine (Ultrasonix Medical Corp. Burnaby, BC, Canada) was used here in the experiments. The computer has an Intel Pentium 4 CPU running at 3.2 GHz and 1 GB of RAM. The center transmission frequency of the transducer was set to 5 MHz, the spatial sampling frequency of the RF signals was 40 MHz, and the RF data capturing depth was set to 5 cm. When the scan plane was set to contain 64 RF scan lines, the frame rate of the RF data was 28 Hz. When the scan plane was set to contain 128 RF scan lines, the RF data frame rate was 14 Hz.  Figure 4.7: A biplane TRUS transducer mounted on a motorized actuator.  To test the proposed algorithms, an agar-agar phantom was built for the experiments, as shown in Fig. 4.7 (b). The phantom consisted of a cubic substrate and a cylindrical hard inclusion at the known position. The substrate was made from water, 2% agar (by mass) and 2% cellulose (A-6924 agar, S-5504 cellulose, Sigma-Aldrich Inc.) as scattering particles. The hard inclusion was made from water, 3% agar and 2% cellulose to create a difference in stiffness between the inclusion and the substrate. The inclusion is approximately as twice stiff as the substrate, as the stiffness of agar is roughly proportional to the squared percentage [47]. The phantom was compressed from the top by a linear probe mounted on a motion stage driven by bandpass filtered white noise between 0.5 and 5.0 Hz. To tune and characterize the TRUS VE system, a prostate phantom was constructed. The phantom consists of a rubber balloon filled with 14% gelatin gel and 2% cellulose (by mass) (G9382B gelatin, S-5504 cellulose, Sigma-Aldrich Inc.). The balloon is embedded in a 10 × 10 × 10 cm3 substrate of 10% gelatin gel and 2% cellulose. The gel in the balloon is approximately twice as stiff as the substrate and has similar echogenic properties. A cylindrical hole in the phantom mimics the rectum and permits insertion of the endorectal ultrasound probe. The prostate phantom was compressed by a biplane TRUS transducer mounted on a motorized actuator (as shown in Fig. 4.7) driven by bandpass filtered white noise between 0.5 and 4.5 Hz. In the experiments, each RF scan line was divided into 80 overlapping segments with length of 71  approximately 1.0 mm. The PSDs and/or CSDs of the displacement estimates of these segments are updated every other RF frame, which means that the frame rate of VE (TF strain) images or power strain images is half of the RF data frame rate. PSD and/or CSD estimation is performed using moving windows with a time sequence of 32 displacement estimates (also the RF frames), and using accumulative windows with a forgetting factor of 0.95. The window length was set to 8 frames, and the window overlap was set to 75% (6 frames). A 16-point FFT was used for computing the PSDs and/or CSDs. Therefore, zero-padding is needed to the data window prior to the FFT. The length for strain LSE was set to 6 pixels. A 2D median filter with a 3 × 3 pixel kernel was applied to the result of VE or power strain images.  4.5  Results  Fig. 4.8 (a)-(c) show the B-mode ultrasound image, VE image and power strain image of the hard inclusion phantom. Fig. 4.8 (d) depicts the center vertical A-lines extracted from the VE image and from the power strain image, respectively. Fig. 4.8 (e) shows the PSD curves at different depths along the center vertical A-line of the power strain image. The contrast and contrast-tonoise ratio (CNR) between the hard inclusion and the background in the VE and power strain images are presented in Table 4.1. The hard inclusion was determined by selecting the pixels that are around the cylinder center and are within the cross-section area of the cylinder. The hard inclusion is invisible in the B-mode image, while it can be clearly seen in the power strain and VE images with high contrast and similar image quality. B-mode and VE images and TF/CF curves of the prostate phantom described above are also shown in Fig. 4.9 (a) - (f).  Table 4.1: Contrast and CNR between the hard inclusion and the background in the VE and power strain images (Fig. 4.8 (b) and (c)) of the hard inclusion phantom.  Contrast CNR  VE image 2.56 9.59  Power strain image 2.59 10.61  Table 4.2: Computational time per frame for power strain imaging and VE using accumulative windows in the phantom experiments. The total time includes the TDPE, PSD and/or TF estimation, and LSE of strain. Number of scan lines 64 128  PSD estimation time (ms) 3.2 6.5  Total time of power strain image (ms) 16 33  TF estimation time (ms) 12 28  Total time of VE image (ms) 25 55  For the computational cost, when the scan plane contains 128 RF scan lines, it takes approximately 59 ms to compute a VE image using TF estimation with moving windows, and it takes 72  (a)  (b)  (c)  1  3 Power strain image VE image 2.5  0.6  2  1/K  PSDs  0.8  0.4  0.2  0  PSD at 4 cm PSD at 3 cm PSD at 2 cm PSD at 1 cm  1.5  1  10  20  30  Depth (mm) (d)  40  50  0.5 0  2  4  6  8  10  12  14  Freq. (Hz) (e)  Figure 4.8: (a) B-mode ultrasound image, (b) VE (TF strain) image and (c) power strain image of the hard inclusion phantom, (d) center vertical A-lines extracted from (b) and (c), and (e) PSD curves at different depths along the center vertical A-line of (c). Color coding of the power strain and VE images is HOT as defined in Matlab (please see the index bar).  about 55 ms to compute a VE image using the TF estimation with accumulative windows. However, when the batch TF estimation is used with a sequence of 32 tissue displacement estimates, it takes nearly 200 ms to compute a VE image. Moreover, Table 4.2 presents the computational time per frame for power strain imaging and VE in the phantom experiments. As one can see, power strain imaging is faster than VE. The computational times are nearly proportional to the number of RF A-lines. The total computing time of the power strain imaging and VE imaging per frame is less than the RF frame period (1/28≈35.7 ms and 1/14≈71.4 ms), which means that the display rate in power strain imaging and VE imaging is limited by the ability of the machine to produce the RF data frames, and not by the computing time needed by the method.  4.6  Discussion  TF and CF estimation with moving windows and with accumulative windows are two fast algorithms for the implementation of real-time ultrasound vibro-elastography. Based on the experimental results, one can see the computational time for TF and CF estimation has been significantly reduced by using moving windows or accumulative windows. The computational load of TF/CF 73  (a)  (b) 2.5  2  1/K  1.5  1  0.5  0  20  40  60  80  Depth  (c)  (d) FR at 4 cm FR at 3 cm FR at 2 cm FR at 1 cm  1  Coherence Function  Frequency Response  1  0.9  0.8  0.7  0.6 0  1  2  3  4  Frequency (Hz) (e)  5  6  7  0.95  0.9  0  CF at 4 cm CF at 3 cm CF at 2 cm CF at 1 cm  1  2  3  4  5  6  7  Frequency (Hz) (f)  Figure 4.9: (a) B-mode image and (b) VE image of the prostate phantom in a transversal scan plan, (c) VE image of the prostate phantom in a sagittal scan plane and (d) relative compliance along the marked VE A-line showing a change of a factor of 2, (e) TFs and (f) corresponding CFs along the marked A-line are displayed. The motion applied is 4 mm peak-to-peak, band-pass filtered (0.5-10 Hz) white noise. Transfer functions are averaged over 10 seconds, and are computed with reference to the motion in the focal area of the firing crystal at a depth of 25 mm. Color coding is HOT as defined in Matlab.  74  estimation mainly focuses on the calculation of FFT-based PSDs and CSDs of tissue motion estimates. Assuming that the number of tissue elements in the scan region is constant, it depends on the length of data windows and the number of FFT points used. The larger the window length or the number of FFT points is, the longer the computational time is. The calculation of TFs and CFs contributes less to the computational cost. Moreover, the computational load is independent of the number of data windows, because the sums of PSDs and CSDs are used and updated every time. This computational speed-up does not involve any loss in quality. If the length of the time sequence of displacement estimates is relatively short, TF estimation with moving windows is better suited for following fast changes in tissue stiffness, because only the PSDs and CSDs of recent data windows are utilized. TF estimation with accumulative windows produces a result that is more accurate and more stable as time goes by, because more historical data are involved. Typically, we choose the forgetting factor λ between 0.90 and 0.98. A small forgetting factor would cause the result image to be inaccurate and unstable, while a large forgetting factor (very close to 1.0) is not able to follow real-time changes in tissue stiffness or changes in probe positions. Furthermore, TF estimation with accumulative windows can be implemented more efficiently and saves more memory compared to TF estimation with moving windows. Thus, TF estimation with accumulative windows is a relatively more flexible choice. The VE and/or power strain images as shown in Figs. 4.8 and 4.9, depict good delineation of the hard inclusion in the phantom by comparison to the B-mode ultrasound image. Furthermore, the results also show high consistency in the computed PSDs and TFs over the frequency range of the excitation, and high coherence over the frequency range of excitation. For good estimation of the TFs and CFs, there must be sufficient frequency content in the transducer motion which provides the tissue compression [114]. In the experiments, the motion estimates of segments in the focal area in the middle of the scan plane were used as a reference. While this motion may not be as “white” as the transducer motion which is under computer control, the transducer location to the RF data acquisition has not been synchronized yet. Such synchronization would be able to improve accuracy and reliability of TF/CF estimation. The power strain image has almost the same image quality and contrast as the VE image, as shown in Table 4.1. The CNR of the power strain image is a little higher than that of the VE image, because VE computation involves a division by the PSD estimates of the reference. However, the PSD estimates of the reference may not be very accurate due to the errors in the tissue motion estimates for the reference segments. Moreover, as the CSD computation is not needed, power strain imaging further reduces the computational cost and saves more memory space compared to the original VE. If it is difficult to accurately estimate the displacements of the reference such as the excitation, the power strain image would be more reliable. As VE provides more complete dynamic tissue description in terms of TFs and CFs, elastic properties of tissue such as stiffness and viscosity could be identified from the TFs, which enables a wide choice of parameters for display. Power strain imaging cannot provide this capability. If the displacements of the reference can be accurately obtained and more elastic properties of tissue are expected, VE is definitely the optimal choice. If only high-quality strain images are expected and a higher frame rate is preferred, power strain imaging is a better option. Additionally, as opposed to conventional free-hand strain 75  imaging techniques, the VE and power strain imaging methods incorporate information about the excitation frequencies. Furthermore, the PSD and/or TF estimation is also equivalent to a temporal averaging procedure, which increases the stability of sequential strain images. Thus, the result would be more reliable, more stable and more robust than conventional strain images. Phantom experiments demonstrate that the VE and power strain imaging techniques produce stable and operator-independent strain images with high signal-to-noise ratio in real time. The experimental results also indicate that the phantom material behaves linearly, as can be seen from the coherence functions. The power strain imaging method can be considered as a simplified and fast version of VE. They may prove useful in clinical applications for prostate brachytherapy, e.g. the detection of prostate contours (as the prostate is a firm organ with respect to its surrounding tissues).  76  Chapter 5  Vibro-Elastography Imaging of the Prostate Region In this chapter, we report our contribution to the 3D vibro-elastography (VE) system that has been designed and built in our research group for imaging the prostate region in patient studies. Patient data have been collected during the volume study (typically 1-2 months in advance) and in the operating room prior to the brachytherapy, at the British Columbia Cancer Agency, Vancouver Centre.  5.1  VE System Design  5.1.1  Ultrasound Machine  Initial data collection was started with an Ultrasonix 500 RP PC-based ultrasound machine, which is the same machine used in Chapters 2 and 4. The maximal depth of RF data capture is 50 mm due to the hardware limitations of that ultrasound machine. Currently the VE imaging system utilizes a Sonix RP PC-based ultrasound machine, as shown in Fig. 5.1, with a dual plane (linear and convex) broadband 5-9 MHz endorectal TRUS transducer (Ultrasonix Medical Corp. Burnaby, BC, Canada). The new Sonix machine does not have the RF data capture limitation of the RP system. The computer has an Intel Pentium 4 CPU running at 3.2 GHz and 3.0 GB of RAM. Both 8-bit B-mode images and 16-bit digitized RF signals can be acquired simultaneously from the machine. The real-time data acquisition program presented in Section 4.3.2 is installed and running on the Sonix machine. The biplane TRUS probe provides two scan planes: the transversal scan plane (a sector image normal to the long axis of the probe) and the sagittal scan plane (a rectangular image containing the longitudinal axis of the probe). The length of the linear crystal array of the transducer is 55 mm. The sector angle of the convex crystal array is 150◦ .  5.1.2  Actuation System  The endorectal transducer is mounted on an actuation system, which is used to apply compression waves to the rectal wall. Two generations of actuation systems have been designed by Salcudean et al. The control program of the VE actuation system is also installed and running on the Sonix machine. The details of the mechanical designs are out of scope in the thesis.  77  Figure 5.1: Photo of the Sonix RP PC-based ultrasound machine (Ultrasonix Medical Corp. Burnaby, BC, Canada).  First-generation Actuation System The first-generation actuation stage for prostate brachytherapy was presented in [114]. Fig. 5.2 shows its schematic drawing, in which the endorectal TRUS probe is mounted in a housing and is vibrated along its sagittal or axial plane by two motor-driven four-bar linkages, along a direction set by the rotational stepper, with its depth adjusted by the translational stepper [114]. Volume images can be acquired using the transverse or the sagittal imaging planes. The actuation stage can be mounted on a standard brachytherapy stand for coarse positioning. Second-generation Actuation System The second-generation actuation system was built in 2008. As illustrated in Fig. 5.3, the endorectal TRUS probe is mounted into the cradle of a CIVCO EX II stepper (CIVCO Medical Solutions, Kalona, IA, USA). With the stepper and cradle, the TRUS transducer can be manually translated along the longitudinal axis of the probe. The cradle has two circular guides at the proximal and distal ends of the TRUS probe. These fit into the metal rollers of the EX II Stepper and allow the cradle to rotate the TRUS transducer about its long axis. The motor and cam mounted on the cradle are used to impart radial vibration to the TRUS probe for VE imaging in the sagittal scan planes. A gear mounted on the cradle engages the shaft of an encoder which allows for accurate 78  Figure 5.2: Schematic drawing of the first-generation VE actuation system. The endorectal TRUS probe is mounted in a housing and is vibrated along its axial plane by two motor-driven four-bar linkages, along a direction set by the rotational stepper, with its depth adjusted by the translational stepper.  measuring of the para-sagittal angle θ of the transducer, when the cradle is manually rotated. In our recent set-up, the original EX II encoder was replaced by a motor with a gear head and its own encoder, allowing for automatic rotation of the transducer through an angle range. The vibration and rotation motors are controlled separately. The rotation angle range can be adjusted via the user interface of the actuator control program. The maximum angle range is between -48◦ and 50◦ , because of mechanical interference. The automatic rotation can be performed in two different modes: (a) Continuous mode: the transducer is smoothly rotated from one end to the other without pauses. The rotation time is 40 seconds; (b) Incremental mode: the transducer is rotated quickly and paused for a short period at certain angles, e.g. pausing 1-5 seconds in 5◦ increments. Typically, the incremental mode is used for pre-operative VE imaging of the prostate region, and the continuous mode is used for post-operative reflected power imaging of the implanted brachytherapy seeds (please refer to Chapter 3). The rotation angle is saved in the header of each RF or B-mode data frame, with a resolution of 0.1◦ .  5.2 5.2.1  3D VE Imaging 2D Signal/Image Processing  When the TRUS transducer is rotated and vibrated automatically at the same time, VE or TF strain images have to be computed based on the angle information with the RF frame headers. When the absolute difference between rotation angles are less than 0.2◦ , the RF data frames are considered to be in the same sagittal scan plane because of physical proximity, angle measuring 79  (a)  (b)  Figure 5.3: (a) Schematic drawing and (b) photo of the second-generation actuation system. The motor and cam mounted on the transducer cradle are used to impart radial vibration to the TRUS probe for VE imaging in the sagittal scan planes.  80  errors and ultrasound beam width. When the probe is rotated in continuous mode (which was not used in our experiments) and when the absolute difference between rotation angles are equal to or greater than 0.5◦ , VE images are considered to be in different para-sagittal scanning planes. Furthermore, when the transducer is rotated in the incremental mode, the first 5 - 10 frames of RF data are discarded in order to skip the transition period from one angle to the next. The remaining RF data are utilized for tissue motion and/or TF estimation. In each sagittal scan plane, tissue displacements due to the probe vibration are determined by the TDPE method [151]. With the time sequences of tissue motion estimates, TFs and CFs between different spatial locations in tissue are estimated by the algorithm presented in Chapter 4, using accumulative windows. The PSD averages of the displacement estimates of the tissue elements in the middle row of the VE imaging region are used as the reference for TF and CF estimation. In each sagittal scan plane, the overall averages of the normalized correlation coefficients and of the CFs are monitored for the reliability of tissue motion and TF estimation, respectively. If the average of correlation coefficients is relatively low, the tissue motion estimates are not reliable and thus the VE image in that sagittal scan plane is discarded. If the average of the CFs is relatively low, the TF estimation results are not reliable and thus the conventional strain image or power strain image will be used instead for that sagittal scan plane.  5.2.2  3D VE Image Reconstruction  The 3D VE image is reconstructed from a series of 2D sagittal VE or strain images computed in Section 5.2.1, when the transducer is rotated and vibrated automatically. When the incremental mode is used, the 2D sagittal VE images are evenly distributed in constant angular increments. Similar to the 3D reflected power image described in Chapter 3, the destination-oriented approach is utilized here for image reconstruction.  5.3  Experiments  Patient data have been collected during the volume study (typically 1 - 2 months before the brachytherapy) and in the operating room (OR) prior to the implant, at the British Columbia Cancer Agency, Vancouver Centre. So far there are 15 patients imaged by the VE imaging system, including the 3 patients that had been imaged in [114], 9 cases recently done during the volume study, and 3 cases recently done in the OR pre-operatively. The Ultrasonix 500 RP machine was used for the first three patients, and the Sonix RP machine was used for the others. Patients 4 - 15 also had the prostate MRI scanned at the UBC hospital, 1 - 2 months before the brachytherapy. The MRI images will be used as a reference for prostate boundary segmentation from 3D VE images. Please refer to Appendices B and C for the subject information and consent forms for the pre-operatively VE imaging during the volume study and in the OR, respectively. Table 5.1 summarizes and compares the data acquisition (DAQ) protocols used for VE imaging of the prostate region during the volume study and in the OR prior to the implant. Please refer to Appendix D for the detailed DAQ protocols. The excitation amplitude (peak to peak) is less than 1% of the strain level with respect to the imaging depth. Because of the limited data collection 81  time and limited space in the OR, the TRUS transducer is automatically rotated and the DAQ duration is much less than that during the volume study. Also in the OR, the imaging depth is decreased if the patient has a small-size prostate, to obtain a higher data frame rate. The excitation frequency range is also increased to contain more periods of the excitation frequency components. Moreover, Table 5.2 presents the parameter values used for tissue motion estimation and VE image computation.  Table 5.1: Summary of the data acquisition (DAQ) protocols used for VE imaging of the prostate region during the volume study and in the OR prior to the implant.  Probe rotation Angle range Angle increment DAQ duration Center frequency RF line count Imaging depth Data frame rate Excitation amplitude Excitation freq range 1 Depending  Volume study Manually -45◦ - 45◦ 5◦ 10 sec 5 MHz 128 60 mm 20 Hz 0.4 - 0.5 mm 0.5 - 4.5 Hz  OR study Automatically -45◦ - 45◦ 1◦ 1 sec 5 MHz 128 50 - 60 mm1 36 - 42 Hz 0.4 - 0.5 mm 2 - 10 Hz  on the size of patients’ prostate  Table 5.2: Parameters used for VE image computation offline. Computed VE image depth Number of RF segments Length of RF segments Spacing between segment centers Length of displacement window Number of FFT points Forgetting factor 1 Depending 2 The  5.4  50 - 60 mm1 100 - 120 1.0 mm 0.5 mm2 8 64 0.95  on the size of patients’ prostate  overlap between RF segments is 50%  Results  Figs. 5.4 (a) - (i) depict the B-mode TRUS images and corresponding VE images of the prostate region of patient 3 in three sagittal scan planes at 0◦ , 10◦ and 20◦ , and the relative compliance curves, four TF curves and the corresponding CF curves along the A-line in Fig. 5.4 (c), at depths of 1, 2, 3 and 4 cm. 82  (a)  (b)  (d)  (e)  (f)  1  2  1/K  1.5  1  0.5  20  40  Depth  (g)  60  80  1  FR at 4 cm FR at 3 cm FR at 2 cm FR at 1 cm  Coherence Function  Frequency Response  2.5  0  (c)  0.9  0.8  0.7  0.6 0  1  2  3  4  Frequency (Hz)  (h)  5  6  7  0.9  0.8 CF at 4 cm CF at 3 cm CF at 2 cm CF at 1 cm 0.7  0  1  2  3  4  5  6  7  Frequency (Hz)  (i)  Figure 5.4: (a) - (c) B-mode TRUS images and (d) - (f) VE images of the prostate region of patient 3 in three sagittal planes at 0◦ , 10◦ and 20◦ , (g) the relative compliance curve, (h) four TF curves and (i) the corresponding CF curves along the A-line in (c), at depths of 1, 2, 3 and 4 cm. Color coding is HOT as defined in Matlab (please see the index bar).  Figs. 5.5 - 5.16 show the B-mode TRUS images, computed VE images, normalized correlation coefficient images with tissue motion estimation, and CF images with TF estimation of the prostate region of patients 4 - 15, respectively, in the sagittal scan planes at different rotation angles. The rotation angles are symmetrical about the vertical (0◦ ). According to the prostate brachytherapy convention, the endorectal transducer is on the bottom of the B-mode and VE images. The left side of the images is close to the apex of the prostate, and the right side is close to the base of the prostate and the bladder (probe inserted from the left). The bladder is on the top-right corner of the images. Prior to the implant, the bladder was injected with certain fluoroscopic contrast dye that is transparent in ultrasound via a catheter in the urethra. For patients 4, 6 and 8, the prostate gland was not entirely visualized because of the inadvertent positioning of the TRUS transducer. For the VE, normalized correlation coefficient and CF images of patients 4 - 15, color coding is 83  JET as defined in Matlab. In the VE images, red color refers to “soft” and blue color corresponds to “stiff”. Additionally, Fig. 5.17 presents the average normalized correlation coefficient and CF curves of patients 4 - 15, at rotation angles between -45◦ and 45◦ in 5◦ increments. The bladder area was cropped manually from the normalized correlation coefficient and CF images.  5.5  Discussion  The VE images obtained not only show good delineation of the prostate gland in most patients by comparison to conventional B-mode TRUS images, but also provide more repeatable details within the prostate region. For instance, the urethra, which is softer than the prostate, can be detected from the VE images of patients 6, 14 and 15 at around 0◦ in Fig. 5.7 (j), Fig. 5.15 (j) and Fig. 5.16 (j) (urethra pointed by a white arrow), while it cannot be found from the corresponding B-mode images. For patient 6, the pubic arch surface can be seen from the VE images (on the top half) as well as from the B-mode images in Fig. 5.7. The calcification inside the prostate gland in patients 9 and 10 can be identified from the VE images in Fig. 5.10 (f) and in Figs. 5.11 (f) and (n) (calcification pointed by a white arrow), as well as from the B-mode images. The calcification and pubic arch appear as dark blue areas in the VE images and as bright areas in the B-mode images. For patient 13, the catheter which was inadvertently left in the urethra during VE imaging is shown as dark blue areas in the VE image in Fig. 5.14 (n) (catheter pointed by a white arrow), as it is much stiffer than surrounding tissue. The prostate delineation and repeatable details in the VE images have been confirmed by a radiation oncologist from the BC Cancer Agency, Vancouver Centre. Thus, prostate segmentation could be improved by the use of VE images. Quantitative evaluation of the resulting VE images of the prostate region was performed by Mahdavi and Moradi et al. in [82]. From Figs. 5.5 - 5.16, one can see that the normalized correlation coefficients associated with displacement estimation are very high in over 90% of the correlation images, except for the bladder region on the top-right corner, the region above the pubic arch as shown in Fig. 5.7, and the calcification region in the prostate as shown in Fig. 5.10 (g) and Figs. 5.11 (g) and (o). This result indicates that the correlation coefficients drop in the region of hard objects during tissue motion estimation. Furthermore, Fig. 5.17 shows that the average normalized correlation coefficients for patients 4 - 15 are above 0.97 at different rotation angles, excluding the bladder area. This indicates that the VE techniques can produce accurate and reliable tissue displacement estimation and resulting VE images. The results in Fig. 5.4 show surprisingly high coherence over the frequencies of excitation, over 0.9 in 80% of TFs shown in the images. While the excitation frequency is quite low, this indicates good confidence in the measurement. In our current work, the PSD averages of the motion estimates of the tissue element along the middle row in the sagittal plane are utilized as a reference for TF and CF estimation. However, the tissue dynamics in the prostate region may not be as “linear” as the tissue model shown in Figs. 1.6. The contact between the TRUS transducer and the rectal wall may not be always good at different angular positions. Also the 84  reference motion may not be as “white” as the transducer motion that is under computer control. Nevertheless, we have not accurately synchronized yet the TRUS transducer position to the RF data acquisition. Such synchronization may improve the accuracy and reliability of TF and CF estimation. Fig. 5.17 shows that the average CFs for patients 4 - 15 are above 0.83 at most rotation angles, excluding the bladder area. Typically, the closer to the vertical or 0◦ , the higher the CF, because the contact between the endorectal probe and the rectal wall may not be good when the transducer is rotated at relatively large angles. Only relative compliance or stiffness can be provided by the VE imaging system presented, as the estimation of absolute values requires the measurement of transducer forces. However, the contrast between the prostate and surrounding tissue appears higher in the relative compliance or stiffness VE images than in the conventional B-mode images. The display of other TF parameters is possible. For example, correlating the changes in TF phase with changes in viscosity may also provide image contrast [32, 114]. Unfortunately, we were not able to collect RF data at sufficiently high frame rate to resolve phase lag, because of frame rate limitations. As mentioned in Chapter 3, there is a signal saturation problem with the Sonix RP machine. The saturation causes sudden phase changes in RF signals when the received RF signal magnitude is very high. The sudden phase change will twist the wave form of RF signals, reduce the correlation coefficients with tissue motion estimation, and cause errors or even false positives in the displacement and TF estimates. We are expecting Ultrasonix Medical Corp. to resolve this problem. Otherwise, we may have to use other ultrasound machines that have a data acquisition interface for RF signals. A B-K ultrasound machine (B-K Medical ApS, Herlev, Denmark) is being used for prostate brachytherapy at the BC Cancer Agency, Vancouver Center. The signal saturation issue does not happen with the B-K machine. The image quality of the B-K machine is also higher than that of the Sonix RP machine. Currently the real-time access of RF signals is not provided by the B-K machine. A frame grabber with a research package will be available in the future for the real-time acquisition of RF data. Overall, the VE imaging technique can be used for prostate segmentation in patients, and the reflected power imaging method (please refer to Chapter 2) has potential in seed detection. Localizing the seed distribution with respect to the prostate boundary would produce the dosimetry for prostate brachytherapy. The VE images and reflected power images can be computed from the same RF data set. When the division of segments in every RF data frame is the same for both VE and reflected power imaging, the VE and reflected power images have the same spatial location and resolution. Thus image registration is not required between the two kinds of images. Seed detection results from a 3D reflected power image (please refer to Chapter 3) and prostate segmentation results from a 3D VE image can be easily combined together to compute the delivered dose. Additionally, image registration between the VE or reflected power images and the B-mode images can be also implemented easily by interpolation as the B-mode image is derived from the same RF data but may have different spatial resolution.  85  (a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)  (i)  (j)  (k)  (l)  (m)  (n)  (o)  (p)  (q)  (r)  (s)  (t)  Figure 5.5: (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 4, in the sagittal planes at -30◦ , -15◦ , 0◦ , 15◦ and 30◦ , from the top down. Color coding is JET as defined in Matlab (please see the index bar). 86  (a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)  (i)  (j)  (k)  (l)  (m)  (n)  (o)  (p)  (q)  (r)  (s)  (t)  Figure 5.6: (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 5, in the sagittal planes at -20◦ , -10◦ , 0◦ , 10◦ and 20◦ , from the top down. 87  (a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)  (i)  (j)  (k)  (l)  (m)  (n)  (o)  (p)  (q)  (r)  (s)  (t)  Figure 5.7: (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 6, in the sagittal planes at -10◦ , -5◦ , 0◦ , 5◦ and 10◦ , from the top down. The white arrow in (j) points to the urethra. 88  (a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)  (i)  (j)  (k)  (l)  (m)  (n)  (o)  (p)  (q)  (r)  (s)  (t)  Figure 5.8: (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 7, in the sagittal planes at -20◦ , -10◦ , 0◦ , 10◦ and 20◦ , from the top down. 89  (a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)  (i)  (j)  (k)  (l)  (m)  (n)  (o)  (p)  (q)  (r)  (s)  (t)  Figure 5.9: (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 8, in the sagittal planes at -20◦ , -10◦ , 0◦ , 10◦ and 20◦ , from the top down. 90  (a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)  (i)  (j)  (k)  (l)  (m)  (n)  (o)  (p)  (q)  (r)  (s)  (t)  Figure 5.10: (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 9, in the sagittal planes at -20◦ , -10◦ , 0◦ , 10◦ and 20◦ , from the top down. The white arrow in (f) points to the calcification. 91  (a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)  (i)  (j)  (k)  (l)  (m)  (n)  (o)  (p)  (q)  (r)  (s)  (t)  Figure 5.11: (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 10, in the sagittal planes at -20◦ , -10◦ , 0◦ , 10◦ and 20◦ , from the top down. The white arrows in (f) and (n) point to the calcification. 92  (a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)  (i)  (j)  (k)  (l)  (m)  (n)  (o)  (p)  (q)  (r)  (s)  (t)  Figure 5.12: (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 11, in the sagittal planes at -10◦ , -5◦ , 0◦ , 5◦ and 10◦ , from the top down. 93  (a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)  (i)  (j)  (k)  (l)  (m)  (n)  (o)  (p)  (q)  (r)  (s)  (t)  Figure 5.13: (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 12, in the sagittal planes at -10◦ , -5◦ , 0◦ , 5◦ and 10◦ , from the top down. 94  (a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)  (i)  (j)  (k)  (l)  (m)  (n)  (o)  (p)  (q)  (r)  (s)  (t)  Figure 5.14: (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 13, in the sagittal planes at -20◦ , -10◦ , 0◦ , 10◦ and 20◦ , from the top down. The white arrow in (n) points to the catheter in the urethra. 95  (a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)  (i)  (j)  (k)  (l)  (m)  (n)  (o)  (p)  (q)  (r)  (s)  (t)  Figure 5.15: (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 14, in the sagittal planes at -30◦ , -15◦ , 0◦ , 15◦ and 30◦ , from the top down. The white arrow in (j) points to the urethra. 96  (a)  (b)  (c)  (d)  (e)  (f)  (g)  (h)  (i)  (j)  (k)  (l)  (m)  (n)  (o)  (p)  (q)  (r)  (s)  (t)  Figure 5.16: (From left to right) B-mode images, VE images, normalized correlation coefficient images and CF images of the prostate region of patient 15, in the sagittal planes at -30◦ , -15◦ , 0◦ , 15◦ and 30◦ , from the top down. The white arrow in (j) points to the urethra. 97  1  0.95 Normalized correlation Coherence function 0.9  0.85  0.8  0.75  −40  −30  −20  −10  0 10 Angle (deg)  20  30  40  Figure 5.17: Average normalized correlation coefficient and average CF curves of patients 4 - 15, excluding the bladder area, at rotation angles between -45◦ and 45◦ in 5◦ increments.  98  Chapter 6  Conclusions In Chapters 2 - 5, we proposed the signal processing methods for enhancing the visibility of brachytherapy seeds, the approach for seed detection in 3D TRUS, the computational algorithms and application programs for real-time vibro-elastography, and the vibro-elastography imaging system for the prostate region. The experimental results and discussion were presented with individual chapters. This chapter concludes with an overview of the research project and a summary of the contributions. The main focus of this research project has been to develop a potential solution towards ultrasound-based dosimetry that is intra-operative, reliable and user-independent for prostate brachytherapy. Such a system would bring many benefits for both patients and healthcare practitioners, e.g. reducing the examination cost and exposure to radiation. While the overall objective is clinically driven, the solution developed during this research addresses it in two parts: seed detection and prostate vibro-elastography imaging, with phantoms and/or patients. Seed detection using ultrasound still remains a challenge in patient studies. Even if only part of the implanted seeds can be located by the proposed seed detection method, the detected seeds together with the detected needle insertion tracks can be used for registration to the fluoroscopic seed detection result, in which fiducial markers are not needed.  6.1  Contributions  The current research has made substantial contributions towards the main goal of the thesis: ultrasound-based, intra-operative dose feedback for prostate brachytherapy. To this end, a potential ultrasound-based intra-operative dosimetry solution has been proposed, and a functional prototype imaging system that makes it possible to test the proposed solution in the lab or clinically has been developed (although the extent of the clinical tests were limited because the Sonix RP ultrasound machine that we used has not been approved for clinical use in prostate brachytherapy, and were limited by a saturation problem in the Sonix RP embedded code that could not be overcome in the timeframe of this thesis). To this end: • We proposed a new ultrasound imaging method, named reflected power imaging, that enhanced the visibility and imaging of brachytherapy seeds. Experimental results have shown that the contrast and contrast-to-noise ratio between the implanted seeds and the background is much higher in the reflected power images than in the conventional B-mode images. The algorithm is not mathematically complex and thus it is suitable for real-time uses. 99  • We proposed a novel solution for automatic detection of brachytherapy seeds in 3D reflected power images, instead of conventional B-mode images, and using a priori knowledge of the needle insertion entries and seed placements along each needle in a pre-operative treatment plan. The computation of a “difference” image before and after the implant is not needed in the proposed method. With the proposed solution, needle insertion tracks (not needles) were detected by using the Hough transform in local 3D spaces, with a priori knowledge of pre-planned needle insertion entries. This is the first time that needle insertion tracks are detected in ultrasound instead of needles. The needle insertion track detection result was used to validate and enhance seed localization along these seed tracks. • The proposed seed detection method was validated using a tissue-equivalent ultrasound prostate phantom, with an implant by a radiation oncologist according to a realistic dosimetry plan generated by a medical physicist. This is the first time that seed detection methods were validated against a realistic treatment plan with such a large number of seeds (there are 136 seeds along 26 needles) in it. The good performance of the proposed solution in the phantom study suggests that it has promise in a clinical setting. • We developed the algorithm for seed segmentation in 2D fluoroscopic images, implemented and validated the fluoroscopic tomosynthesis approach to reconstruct the seed distribution in 3D spaces, from the seed segmentation result in 2D fluoroscopic images at five different perspectives. Not only the centroid but also the orientation of brachytherapy seeds were reconstructed from our tomosynthesis program. The fluoroscopic tomosynthesis result can be used to validate ultrasound-based seed detection results intra-operatively, and can be directly used as a backup plan for dose computation, because fluoroscopy is clinically used during prostate brachytherapy as well. • We implemented fast algorithms for the estimation of transfer functions and coherence functions by using Welch’s method, for the implementation of vibro-elastography in real time. A simplified version of vibro-elastography, called power strain imaging, was also proposed to obtain higher frame rates and high-quality strain images. Moreover, two GUI-based Windows application programs have been designed and developed for RF data analysis and for real-time data acquisition, respectively. The RF data analysis program can compute and display a number of 2D image types (reflected power images, strain images, vibro-elastography images, etc.) from the RF data set, and can store the result data and images for further image analysis. The real-time data acquisition program has been successfully used for phantom data collection in the lab and for patient data collection at the British Columbia Cancer Agency, Vancouver Centre. • A vibro-elastography imaging system has been developed, including the robotic actuation system designed by Salcudean et al. The above-mentioned algorithms for the reflected power imaging and vibro-elastography, components and programs have been integrated into the vibro-elastography system. It has been successfully utilized for real-time data collection of 100  the prostate region of patients. Over 15 patients have been scanned during the volume study and in the operating room prior to the implant, at the British Columbia Cancer Agency, Vancouver Centre, and the results are encouraging. Several peer-reviewed publications have arisen from this work, including one journal paper: • X. Wen et al , Detection of Brachytherapy Seeds Using 3D Transrectal Ultrasound, accepted by the IEEE Transactions on Biomedical Engineering, May 2010, which presents the proposed seed detection solution in 3D reflected power images with a priori knowledge in Chapter 3 and the experimental results with a prostate phantom according to a realistic treatment plan. and seven conference and workshop presentations: • X. Wen et al , Imaging of the Prostate with Vibro-elastography: Preliminary Patient Results, the Workshop on Prostate Image Analysis and Computer-Assisted Intervention (in conjunction with MICCAI 2008), Sept. 2008, which presents the preliminary patient results of vibro-elastography imaging of the prostate region. • X. Wen et al , Detection of Brachytherapy Seeds using 3D Ultrasound, the 30th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, Aug. 2008, which presents the old-version seed detection method in 3D TRUS and the detection result of part of the implanted seeds in the CIRS prostate phantom. • X. Wen et al , Ultrasound-only Dosimetry for Prostate Brachytherapy: Preliminary Phantom Results, the SPIE International Symposium on Medical Imaging, Feb. 2007, which presents the seed detection and hard-inclusion imaging results of a manually-made phantom containing 16 dummy seeds. • X. Wen et al , Power strain imaging based on vibro-elastography techniques, the SPIE International Symposium on Medical Imaging, Feb. 2007, which presents the implementation of fast transfer function estimation for vibro-elastography and a simplified version named power strain imaging. • X. Wen et al , Real-time Kalman-filtered strain imaging of the human prostate, the 29th Canadian Medical and Biological Engineering Conference, June 2006, which presents the preliminary patient results of strain imaging of the prostate region. • X. Wen et al , Fast computation for ultrasound vibro-elastography, the 29th Canadian Medical and Biological Engineering Conference, June 2006, which presents the preliminary implementation of fast transfer function estimation for vibro-elastography. • X. Wen et al , Detection of brachytherapy seeds using ultrasound radio frequency signals, the SPIE International Symposium on Medical Imaging, Feb. 2006, which presents the reflected power imaging method for enhancing seed visibility. Additionally, one patent application has been filed from this work: 101  • X. Wen et al , Method and apparatus for processing ultrasound image signals, US Patent Application, in Feb 2007, which presents the reflected power, squared envelope and decompressed B-mode imaging methods.  6.2  Future Work  As with many research projects, areas that require further work always remain. This section includes some recommendations for future work. The Sonix RP ultrasound machine that we are using in the experiments has some limitations. A B-K ultrasound machine (B-K Medical ApS, Herlev, Denmark) is being used for prostate brachytherapy at the British Columbia Cancer Agency, Vancouver Center. A frame grabber with a research package will be available in the future for the real-time acquisition of RF data. We should migrate the reflected power and vibro-elastography imaging system to the B-K ultrasound machine so that intra-operative patient data can be acquired during the implant in the operating room. Reflected power images and vibro-elastography images computed with the RF data from the B-K machine will be tested for seed detection and prostate segmentation. Furthermore, the image quality of the B-K machine is higher than that of the Sonix RP machine, and the signal saturation issue does not happen with the B-K machine. For seed detection in patient studies, typically needle insertion tracks are not straight lines. They can be approximated either by parabolas or by a series of segments with different orientations. As a result, the searching space for each needle entry will not be a simple cylinder and the Hough transform will also be more mathematically complex. Moreover, if the frame grabber and research package for the B-K ultrasound machine are ready or if the Sonix RP machine gets approved for brachytherapy uses, we can collect RF data every time a radiation oncologist performs a few needle entries. In this way, only the implanted seeds by the most recent needle entries need to be detected every time. As the search volume and the number of seeds to be identified are limited, the probability of finding all the implanted seeds increases. Furthermore, registration between the detected needle insertion tracks from 3D ultrasound and from fluoroscopic tomosynthesis will be implemented as a backup solution for clinical seed detection. For vibro-elastography imaging, RF data acquisition needs to be synchronized with TRUS transducer locations in real time. Such synchronization would be able to improve accuracy and reliability of the transfer function and coherence function estimation. Moreover, if possible, higher RF frame rates are expected to obtain high-frequency components of tissue motion, and to enable the display of other transfer function parameters. Prostate segmentation results from vibroelastography images will be quantitatively evaluated against the manual segmentation results from B-mode images and from MRI images by radiation oncologists. Finally, as this research is clinically driven, more patient studies will need to be performed to validate and approve the feasibility and reliability of the proposed solution and system.  102  Bibliography [1] Canadian cancer statistics 2009. Technical report, Canadian Cancer Society, Inc., Toronto, ON, 2009. [2] Cancer facts and figures 2009. Technical report, American Cancer Society, Inc., Atlanta, GA, 2009. [3] R. G. Aarnink, S. D. Pathak, and et al. Edge detection in prostatic ultrasound images using integrated edge maps. Ultrasonics, 36:635–642, 1998. [4] P. Abolmaesumi, S. E. Salcudean, W. H. Zhu, and et al. Image-guided control of a robot for medical ultrasound. IEEE Trans. Robot. Automat., 18(1):11–23, 2002. [5] P. Abolmaesumi and M. R. Sirouspour. An interacting multiple model probabilistic data association filter for cavity boundary extraction from ultrasound images. IEEE Trans. Med. Imag., 23(6):772–784, 2004. [6] P. Abolmaesumi and M. R. Sirouspour. Segmentation of prostate contours from ultrasound images. In Proc. IEEE ICASSP, volume 3, pages 517–520, 2004. [7] S. K. Alam, E. J. Feleppa, A. Kalisz, and et al. Prostate elastography - Preliminary in vivo results. In Proc. SPIE, volume 5750, pages 339–345, 2005. [8] S. K. Alam and J. Ophir. Reduction of signal decorrelation from mechanical compression of tissues by temporal stretching: Applications to elastography. Ultrason. Imaging, 23:95–105, 1997. [9] R. Alterovitz, J. Pouliot, R. Taschereau, I. Hsu, and K. Goldberg. Needle insertion and radioactive seed implantation in human tissues: Simulation and sensitivity analysis. In Proc. IEEE ICRA, volume 2, pages 1793–1799, 2003. [10] R. Alterovitz, J. Pouliot, R. Taschereau, I. Hsu, and K. Goldberg. Sensorless planning for medical needle insertion procedures. In Proc. IEEE/RSJ IROS, volume 4, pages 3337–3343, 2003. [11] V. G. Andreev, A. E. Parnomaryov, and et al. Detection of prostate cancer with optoacoustic tomography: Feasibility and modeling. In Proc. SPIE, volume 4960, pages 45–57, 2003. [12] J. Awad, T. K. Abdel-Galil, M. M. A. Salama, and et al. Prostate’s boundary detection in transrectal ultrasound images using scanning technique. In Proc. IEEE CCECE, volume 2, pages 1199–1202, 2003. [13] S. Badiei. Prostate segmentation in ultrasound images using image warping and ellipsoid fitting. Master’s thesis, Univ. B.C., Vancouver, 2007. [14] S. Badiei, S. E. Salcudean, J. Varah, and W. J. Morris. Prostate segmentation in 2D ultrasound images using image warping and ellipse fitting. In Proc. MICCAI, volume 4191, pages 17–24, 2006. 103  [15] T. Baldeweck, P. Laugier, and et al. Application of autoregressive spectral analysis for ultrasound attenuation estimation: Interest in highly attenuating medium. IEEE Trans. Ultrason., Ferroelect., Freq. Contr., 42(1):99–110, 1995. [16] L. Beaulieu, D. Tubic, J. Pouliot, and et al. Post-implant dosimetry using fusion of ultrasound images with 3D seed coordinates from fluoroscopic images in transperineal interstitial permanent prostate brachytherapy. Int. J. Radiat. Oncol. Biol. Phys., 48(3 suppl.). [17] V. Behar, D. Adam, and et al. Improving motion estimation by accounting for local image distortion. Ultrasonics, 43:57–65, 2004. [18] J. S. Bendat and A. G. Piersol. Random Data: Analysis and Measurement Procedures. Wiley Interscience, Maryland, USA, second edition, 1986. [19] J. Bercoff, S. Chaffai, and et al. In vivo breast tumor detection using transient elastography. Ultrasound Med. Biol., 29:1387–1396, 2003. [20] I. Cespedes and J. Ophir. Reduction of image noise in elastography. Ultrason. Imaging, 15:89–102, 1993. [21] P. Chaturvedi, M. F. Insana, and T. J. Hall. 2-D companding for noise reduction in strain imaging. IEEE Trans. Ultrason., Ferroelect., Freq. Contr., 45:179–191, 1998. [22] C. H. Chen, J. Y. Lee, W. H. Yang, C. M. Chang, and Y. N. Sun. Segmentation and reconstruction of prostate from transrectal ultrasound images. Biomed. Eng., 8:287–292, 1996. [23] D. Christensen. Ultrasonic Bioinstrumentation. John Wiley and Sons, New Jersey, 1988. [24] D. C. Crawford, D. S. Dell, and J. C. Bamber. Compensation for the signal processing characteristics of ultrasound B-mode scanners in adaptive speckle reduction. Ultrasound Med. Biol., 19(6):469–485, 1993. [25] B. J. Davis, R. R. Kinnick, M. Fatemi, and et al. Measurement of the ultrasound backscatter signal from three seed types as a function of incidence angle: Application to permanent prostate brachytherapy. Int. J. Radiat. Oncol. Biol. Phys., 57:1174–1182, 2003. [26] C. de Korte, E. I. Cespedes, and et al. Intravascular elasticity imaging using ultrasound: Feasibility studies in phantoms. Ultrasound Med. Biol., 23:735–746, 1997. [27] E. Dehghan, X. Wen, R. Zahiri-Azar, M. Marchal, and S. E. Salcudean. Needle-tissue interaction modeling using ultrasound-based motion estimation: Phantom study. Computer Aided Surgery, 13(5):265–280, 2008. [28] M. Ding, Z. Wei, D. B. Downey, and A. Fenster. Automated seed localization for introoperative prostate brachytherapy based on 3D line segment patterns. In Proc. SPIE, volume 5744, pages 417–424, 2005. [29] M. Ding, Z. Wei, L. Gardi, D. B. Downey, and A. Fenster. Needle and seed segmentation in intro-operative 3D ultrasound-guided prostate brachytherapy. Ultrasonics, 44:331–336, 2006. [30] V. Dutt and J. F. Greenleaf. Statistics of the log-compressed echo envelope. J. Acoust. Soc. Am., 99(6):3817–3825, 1996. [31] R. Ebrahimi, S. Okazawa, R. Rohling, and S. E. Salcudean. Hand-held steerable needle device. In Proc. MICCAI, volume 2879, pages 223–230, 2003. 104  [32] H. Eskandari, S. E. Salcudean, and R. Rohling. Viscoelastic parameter estimation based on spectral analysis. IEEE Trans. Ultrason., Ferroelect., Freq. Contr., 55(7):1611–1625, 2008. [33] M. Fatemi and J. Greenleaf. Probing the dynamics of tissue at low frequencies with the radiation force of ultrasound. Phys. Med. Biol., 45:1449–1464, 2000. [34] E. J. Feleppa, S. Ramachandran, S. K. Alam, and et al. Novel methods of analyzing radiofrequency echo signals for the purpose of imaging brachytherapy seeds used to treat prostate cancer. In Proc. SPIE, volume 4687, pages 127–138, 2002. [35] A. Fenster, S. Tong, H. N. Cardinal, C. Blake, and D. B. Downey. Three-dimensional ultrasound imaging system for prostate cancer diagnosis and treatment. IEEE Trans. Instrum. Meas., 47(6):1439–1447, 1998. [36] A. Fort, C. Manfredi, and S. Rocchi. Recursive autoregressive spectral maps for ocular pathology detection. Ultrasound Med. Biol., 23(3):391–403, 1997. [37] D. French. Real-time dosimetry for prostate brachytherapy using TRUS and fluoroscopy. Master’s thesis, Univ. B.C., Vancouver, 2004. [38] D. French, J. Morris, M. Keyes, O. Goksel, and S. E. Salcudean. Computing intraoperative dosimetry for prostate brachytherapy using TRUS and fluoroscopy. Acad. Radiol., 12:1262– 1272, 2005. [39] D. French, J. Morris, M. Keyes, and S. E. Salcudean. Real-time dosimetry for prostate brachytherapy using TRUS and fluoroscopy. In Proc. MICCAI, volume 3217, pages 983–991, 2004. [40] L. Gao, K. J. Parker, R. M. Lerner, and et al. Theory and application of sonoelasticity imaging. Int. J. Imag. Sys. Tech., 8:104–109, 1997. [41] B. S. Garra, E. I. Cespedes, J. Ophir, and et al. Elastography of breast lesions: Initial clinical results. Radiology, 202:79–86, 1997. [42] A. Ghanei, H. Soltanian-Zadeh, A. Ratkewicz, and F. F. Yin. A three-dimensional deformable model for segmentation of human prostate from ultrasound images. Med. Phys., 28:2147– 2153, 2001. [43] J. M. Girault, F. Ossant, and et al. Time-varying autoregressive spectral estimation for ultrasound attenuation in tissue characterization. IEEE Trans. Ultrason., Ferroelect., Freq. Contr., 45(3):650–659, 1998. [44] L. Gong, P. S. Cho, B. H. Han, and et al. Ultrasonography and fluoroscopic fusion for prostate brachytherapy dosimetry. Int. J. Radiat. Oncol. Biol. Phys., 54:1322–1330, 2002. [45] L. Gong, S. Pathak, D. Haynor, P. Cho, and K. Kim. Parametric shape modeling using deformable superellipses for prostate segmentation. IEEE Trans. Med. Imag., 23:340–349, 2004. [46] R. C. Gonzales and R. E. Woods. Digital Image Processing. Prentice-Hall, New Jersey, second edition, 2002. [47] T. J. Hall, M. Bilgen, M. F. Insana, and T. A. Krouskop. Phantom materials for elastography. IEEE Trans. Ultrason., Ferroelect., Freq. Contr., 44(6):1355–1365, 1997. 105  [48] B. H. Han, K. Wallner, G. Merrick, and et al. Prostate brachytherapy seed identification on post-implant TRUS images. Med. Phys., 30:898–900, 2003. [49] K. M. Hiltawsky and et al. ultrasound elastography of breast lesions: Clinical results. Ultrasound Med. Biol., 27:1461–1469, 2001. [50] A. C. Hodge, A. Fenster, D. B. Downeyb, and H. M. Ladak. Prostate boundary segmentation from ultrasound images using 2D active shape models: Optimisation and extension to 3D. Comput. Methods Programs Biomed., 84:99–113, 2006. [51] D. R. Holmes, B. J. Davis, C. J. Bruce, and R. A. Robb. 3D visualization, analysis, and treatment of the prostate using trans-urethral ultrasound. Computerized Med. Imag. Graphics, 27:339–349, 2003. [52] D. R. Holmes, B. J. Davis, C. J. Bruce, T. Wilson, and R. A. Robb. Trans-urethral ultrasound imaging of the prostate for applications in prostate brachytherapy: Analysis of phantom and in vivo data. In Proc. SPIE, volume 4319, pages 46–52, 2001. [53] D. R. Holmes and R. A. Robb. Improved automated brachytherapy seed localization in trans-urethral ultrasound data. In Proc. SPIE, volume 5367, pages 353–360, 2004. [54] Hough transform. Wikipedia, the free encyclopedia, http://en.wikipedia.org/. [55] A. Jain, T. Mustufa, Y. Zhou, C. Burdette, G. S. Chirikjian, and G. Fichtinger. FTRAC A robust fluoroscope tracking fiducial. Med. Phys., 32:3185–3198, 2005. [56] A. Jain, Y. Zhou, T. Mustufa, E. C. Burdette, G. S. Chirikjian, and G. Fichtinger. Matching and reconstruction of brachytherapy seeds using the Hungarian algorithm (MARSHAL). Med. Phys., 32:1153–1171, 2005. [57] A. Jain, Y. Zhou, T. Mustufa, E. C. Burdette, G. S. Chirikjian, and G. Fichtinger. Matching and reconstruction of brachytherapy seeds using the Hungarian algorithm (MARSHAL). In Proc. SPIE, volume 5744, pages 810–821, 2005. [58] F. Kallel and J. Ophir. A least-squares strain estimator for elastography. Ultrason. Imaging, 19:195–208, 1997. [59] F. Kallel and J. Ophir. Three-dimensional tissue motion and its effect on image noise in elastography. IEEE Trans. Ultrason., Ferroelect., Freq. Contr., 44:1286–1296, 1997. [60] F. Kallel, R. R. Price, E. Konofagou, and J. Ophir. Elastographic imaging of the normal canine prostate in vitro. Ultrason. Imaging, 21:201–215, 1999. [61] D. Kaplan and Q. Ma. On the statistical characteristics of log-compressed Rayleigh signals: Theoretical formulation and experimental results. J. Acoust. Soc. Am., 95(3):1396–1400, 1994. [62] A. Klauser and et al. Real-time elastography for prostate cancer detection. J. Urology, 171(4 Suppl.):477, 2004. [63] C. Knoll, M. Alcaniz, and et al. Multiresolution segmentation of medical images using shape-restricted snakes. In Proc. SPIE, pages 222–233, 1999. [64] C. Knoll, M. Alcaniz, V. Grau, C. Monserrat, and M. C. Juan. Outlining of the prostate using snakes with shape restrictions based on the wavelet transform. Pattern Recognition, 32:1767–1781, 1999. 106  [65] K. Konig, U. Scheipers, A. Pesavento, and et al. Initial experiences with real-time elastography guided biopsies of the prostate. J. Urology, 174:115–117, 2005. [66] E. Konofagou and J. Ophir. A new elastographic method for estimation and imaging of lateral displacements, lateral strains, corrected axial strains and poisson’s ratios in tissues. Ultrasound Med. Biol., 24:1183–1199, 1998. [67] E. E. Konofagou, J. D’Hooge, and J. Ophir. Myocardial elastography - A feasibility study in vivo. Ultrasound Med. Biol., 28:475–482, 2002. [68] E. E. Konofagou, T. Harrigan, and S. Solomon. Assessment of regional myocardial strain using cardiac elastography: Distinguishing infarcted from non-infarcted myocardium. In Proc. IEEE Ultrasonics Symp, volume 2, pages 1589–1592, 2001. [69] E. E. Konofagou, T. Varghese, J. Ophir, and et al. Power spectral strain estimators in elastography. Ultrasound Med. Biol., 25:1115–1129, 1999. [70] T. A. Krouskop, T. M. Wheeler, F. Kallel, and et al. Elastic moduli of breast and prostate tissues under compression. Ultrason. Imaging, 20(4):260–274, 1998. [71] S. Kruse, J. Smith, A. Lawrence, and et al. Tissue characterization using magnetic resonance elastography: Preliminary results. Phys. Med. Biol., 45:1579–1590, 2000. [72] H. M. Ladak, F. Mao, Y. Q. Wang, D. B. Downey, A. Steinman, and D. A. Fenster. Prostate segmentation from 2D ultrasound images. Med. Phys., 27:1777–1788, 2000. [73] V. Lagerburg and et al. Measurement of prostate rotation during insertion of needles for brachytherapy. Radiothera. Oncol., 77:318–323, 2005. [74] S. Langeland, J. d’Hooge, H. Torp, B. Bijnens, and P. Suetens. Comparison of time-domain displacement estimators for twodimensional RF tracking. Ultrasound Med. Biol., 29:1177– 1186, 2003. [75] R. Lerner, K. Parker, J. Holen, R. Gramiak, and R. Waag. Sonoelasticity in ultrasonic tissue characterization and echographic imaging. Acoust. Imaging, 19:317–327, 1988. [76] S. Levinson, M. Shinagawa, and T. Sato. Sonoelastic determination of human skeletal muscle elasticity. J. Biomech., 28:1145–1154, 1995. [77] L. Ljung. System Identification, Theory for the User. Prentice Hall, New Jersey, 1999. [78] A. Lorenz, C. Haas, and H. Ermert. Segmentation of ultrasonic prostate images using a probabilistic model based on Markov random processes. Ultrason. Imaging, 19:44–45, 1997. [79] A. Lorenz, G. Schapers, H. J. Sommerfeld, and et al. On the use of a modified optical flow algorithm for the correction of axial strain estimates in ultrasonic elastography for medical diagnosis. In Proc. Forum Acusticum, 1999. [80] A. Lorenz, H. J. Sommerfeld, M. Garcia-Schurmann, and et al. A new system for the acquisition of ultrasonic multicompression strain images of the human prostate in vivo. IEEE Trans. Ultrason., Ferroelect., Freq. Contr., 46:1147–1153, 1999. [81] M. A. Lubinski, S. Y. Emelianov, and M. O’Donnell. Speckle tracking methods for ultrasonic elasticity imaging using shorttime correlation. IEEE Trans. Ultrason., Ferroelect., Freq. Contr., 46:82–96, 1999. 107  [82] S. S. Mahdavi, M. Moradi, X. Wen, W. Morris, and S. Salcudean. Vibro-elastography for visualization of the prostate region: Method evaluation. In Proc. MICCAI, volume 5762, pages 339–347, 2009. [83] S. S. Mahdavi, W. Morris, I. Spadinger, N. Chng, O. Goksel, and S. Salcudean. 3D prostate segmentation in ultrasound images based on tapered and deformed ellipsoids. In Proc. MICCAI, volume 5762, pages 960–967, 2009. [84] J. Mamou and E. J. Feleppa. Singular spectrum analysis applied to ultrasonic detection and imaging of brachytherapy seeds. J. Acoust. Soc. Am., 121:1790–1801, 2007. [85] J. Mamou, S. Ramachandran, and E. J. Feleppa. Angle-dependent ultrasonic detection and imaging of brachytherapy seeds using singular spectrum analysis. J. Acoust. Soc. Am., 123(4):2148–2159, 2008. [86] J. Mamou, S. Ramachandran, and E. J. Feleppa. Angle-dependent ultrasonic detection and imaging of two types of brachytherapy seeds using singular spectrum analysis. J. Acoust. Soc. Am., 124(6):EL347–EL352, 2009. [87] S. A. McAleavey. Doppler technique for the detection and localization of modified brachytherapy seeds. In Proc. IEEE Ultrasonics Symp, volume 4687, pages 190–198, 2002. [88] S. A. McAleavey, M. Palmeri, S. Gracewski, and G. Trahey. Ferromagnetic brachytherapy seed motion in soft tissue: Models, measurements and ultrasound detection. In Proc. SPIE, volume 2, pages 1575–1579, 2002. [89] S. A. McAleavey, D. J. Rubens, and K. J. Parker. Doppler ultrasound imaging of magnetically vibrated brachytherapy seeds. IEEE Trans. Biomed. Eng., 50:252–255, 2003. [90] F. G. Mitri, B. J. Davis, J. F. Greenleaf, and M. Fatemi. In vitro comparative study of vibroacoustography versus pulse-echo ultraosund in imaging permanent prostate brachytherapy seeds. Ultrasonics, 49:31–38, 2009. [91] F. G. Mitri, P. Trompette, and J. Y. Chapelon. Using vibro-acoustography to detect brachytherapy metal seeds. In Proc. IEEE Ultrasonics Symp, volume 2, pages 1528–1531, 2003. [92] F. G. Mitri, P. Trompette, and J. Y. Chapelon. Improving the use of vibro-acoustography for brachytherapy metal seed imaging: A feasibility study. IEEE Trans. Med. Imag., 23:1–6, 2004. [93] S. Nag, W. Bice, K. DeWyngaert, and et al. The American Brachytherapy Society recommendations for permanent prostate brachytherapy postimplant dosimetric analysis. Int. J. Radiat. Oncol. Biol. Phys., 46:221–230, 2000. [94] S. Nag, J. P. Ciezki, R. Cormack, and et al. Intraoperative planning and evaluation of permanent prostate brachytherapy: Report of the American Brachytherapy Society. Int. J. Radiat. Oncol. Biol. Phys., 51:1422–1430, 2001. [95] R. K. Nightingale, M. L. Palmeri, R. W. Nightingale, and G. E. Trahey. On the feasibility of remote palpation using acoustic radiation force. J. Acoust. Soc. Am., 110(1):625–634, 2001. [96] J. A. Noble and D. Boukerroui. Ultrasound image segmentation: A survey. IEEE Trans. Med. Imag., 25:987–1010, 2006. 108  [97] M. O’Donnell, A. R. Skovoroda, B. M. Shapo, and S. Y. Emelianov. Internal displacement and strain imaging using ultrasonic speckle tracking. IEEE Trans. Ultrason., Ferroelect., Freq. Contr., 41:314–325, 1994. [98] J. Ophir, I. Cespedes, H. Ponnekanti, Y. Yazdi, and X. Li. Elastography: A method for imaging the elasticity of biological tissues. Ultrason. Imaging, 13:111–134, 1991. [99] A. Oppenheim, R. Schafer, and J. Buck. Discrete-Time Signal Processing. Prentice-Hall, New Jersey, second edition, 1999. [100] P. F. Orio, I. B. Tutar, S. Narayanan, and et al. Intraoperative TRUS-fluoroscopy fusion improves prostate brachytherapy quality. Int. J. Radiat. Oncol. Biol. Phys., 69:302–307, 2007. [101] N. Otsu. A threshold selection method from gray-level histograms. IEEE Trans. Syst., Man, Cybern., 9:62–66, 1979. [102] S. D. Pathak, R. Aarnink, Y. Kim, and et al. Quantitative three-dimensional transrectal ultrasound (TRUS) for prostate imaging. In Proc. SPIE, volume 3335, pages 83–92, 1998. [103] S. D. Pathak, V. Chalana, D. R. Haynor, and Y. Kim. Edge-guided boundary delineation in prostate ultrasound images. IEEE Trans. Med. Imag., 19:1211–1219, 2000. [104] A. Pesavento and A. Lorenz. Real time strain imaging and in-vivo applications in prostate cancer. In Proc. IEEE Ultrasonics Symp., volume 2, pages 1647–1652, 2001. [105] A. Pesavento, A. Lorenz, S. Siebers, and H. Ermert. New real-time strain imaging concepts using diagnostic ultrasound. Phys. Med. Biol., 45:1423–1435, 2000. [106] A. Pesavento, C. Perrey, M. Krueger, and H. Ermert. A time efficient and accurate strain estimation concept for ultrasonic elastography using iterative phase zero estimation. IEEE Trans. Ultrason., Ferroelect., Freq. Contr., 46:1057–1067, 1999. [107] Placement of seeds for prostate brachytherapy - Radiotherapy Photos. http:// radiotherapy.gemzies.com/show/entry 10085/placement of seeds for prosta.html. [108] L. Potters. Permanent prostate brachytherapy in men with clinically localized prostate cancer. Clin. Oncol., 15:301–315, 2003. [109] R. W. Prager, A. H. Gee, G. M. Treece, and et al. Decompression and speckle detection for ultrasound images using the homodyned K-distribution. Pattern Recognition Letters, 24:705–713, 2003. [110] J. S. Prater and W. D. Richard. Segmenting ultrasound images of the prostate using neural networks. Ultrason. Imaging, 14:159–185, 1992. [111] B. R. Prestidge, D. C. Hoak, P. D. Grimm, H. Ragde, W. Cavanagh, and J. C. Blasko. Posttreatment biopsy results following interstitial brachytherapy in early-stage prostate cancer. Int. J. Radiat. Oncol. Biol. Phys., 37:31–39, 1997. [112] J. B. Priestly and D. C. Beyer. Guided brachytherapy for the treatment of confined prostate cancer. Urology, 40:27–32, 1992. [113] M. J. Rivard, B. M. Coursey, L. A. DeWerd, and et al. Update of AAPM task group No. 43 report: A revised AAPM protocol for brachytherapy dose calculations. Med. Phys., 31:633–674, 2004. 109  [114] S. E. Salcudean, D. French, S. Bachmann, and et al. Viscoelasticity modeling of the prostate region using vibro-elastography. In Proc. MICCAI, volume 4190, pages 389–396, 2006. [115] S. E. Salcudean, T. D. Prananta, W. J. Morris, and I. Spadinger. A robotic needle guide for prostate brachytherapy. In Proc. IEEE ICRA, pages 2975–2981, 2008. [116] L. Sandrin, S. Catheline, M. Tanter, X. Hennequin, and M. Fink. Timeresolved pulsed elastography with ultrafast ultrasonic imaging. Ultrason. Imaging, 21:259–272, 1999. [117] F. Shao, K. V. Ling, W. S. Ng, and R. Y. Wu. Prostate boundary detection from ultrasonographic images. J. Ultrasound Med., 22(6):605–623, 2003. [118] D. Shen, Y. Zhan, and C. Davatzikos. Segmentation of prostate boundaries form ultrasound images using statistical shape model. IEEE Trans. Med. Imag., 23:539–551, 2003. [119] V. Singh, L. Mukherjee, J. Xu, K. R. Hoffmann, P. M. Dinu, and M. Podgorsak. Brachytherapy seed localization using geometric and linear programming techniques. IEEE Trans. Med. Imag., 26(9):1291–1304, 2007. [120] R. Sinkus, J. Lorenzen, D. Schrader, and et al. Tissue characterization using magnetic resonance elastography: Preliminary results. Phys. Med. Biol., 45:1649–1664, 2000. [121] A. Skovoroda, S. Emelianov, and M. O’Donnell. Tissue elasticity reconstruction based on ultrasonic displacement and strain images. IEEE Trans. Ultrason., Ferroelect., Freq. Contr., 42:747–765, 1995. [122] R. Souchon, J. Chapelon, and et al. Feasibility of monitoring HIFU prostate cancer therapy using elastography. In Proc. SPIE, volume 4325, pages 392–399, 2001. [123] R. Souchon, O. Rouviere, and et al. Visualisation of HIFU lesions using elastography of the human prostate in vivo: Preliminary results. Ultrasound Med. Biol., 29:1007–1015, 2003. [124] S. Srinivasan and J. Ophir. A zero-crossing strain estimator for elastography. Ultrasound Med. Biol., 29:227–238, 2003. [125] R. G. Stock, N. N. Stone, and C. Iannuzzi. Sexual potency following interactive ultrasoundguided brachytherapy for prostate cancer. Int. J. Radiat. Oncol. Biol. Phys., 35:267–272, 1996. [126] S. H. Stokes, J. D. Real, P. W. Adams, J. C. Clements, S. Wuertzer, and W. Kan. Transperineal ultrasound-guided radioactive seed implantation for organ-confined carcinoma of the prostate. Int. J. Radiat. Oncol. Biol. Phys., 37:337–341, 1997. [127] Y. Su, B. J. Davis, M. G. Herman, and R. A. Rob. Fluoroscopy to ultrasound image registration using implanted seeds as fiducials during permanent prostate brachytherapy. In Proc. SPIE, volume 5367, pages 371–378, 2004. [128] H. Talhami, L. Wilson, and et al. Power spectral strain estimators in elastography. Ultrasound Med. Biol., 20:759–772, 1994. [129] D. A. Todor, G. N. Cohen, H. I. Amols, and M. Zaider. Operator-free, film-based 3D seed reconstruction in brachytherapy. Phys. Med. Biol., 47:2031–2048, 2002. [130] D. A. Todor, M. Zaider, G. N. Cohen, M. F. Worman, and M. J. Zelefsky. Intraoperative dynamic dosimetry for prostate implants. Phys. Med. Biol., 48:1153–1171, 2003. 110  [131] S. Tong, D. B. Downey, H. N. Cardinal, and A. Fenster. A three-dimensional ultrasound prostate imaging system. Ultrasound Med. Biol., 22:735–746, 1996. [132] A. Tornes and M. Eriksen. A new brachytherapy seed design for improved ultrasound visualization. In Proc. IEEE Ultrasonics Symp, volume 2, pages 1278–1283, 2003. [133] D. Tubic, A. Zaccarin, L. Beaulieu, and J. Pouliot. Automated seed detection and threedimensional reconstruction II: Reconstruction of permanent prostate implants using simulated annealing. Phys. Med. Biol., 28(1):2272–2286, 2001. [134] D. Tubic, A. Zaccarin, J. Pouliot, and L. Beaulieu. Automated seed detection and threedimensional reconstruction I: Seed localization from fuoroscopic images or radiographs. Phys. Med. Biol., 28(1):2265–2271, 2001. [135] E. Turgay. Imaging visco-elastic properties of soft tissue with ultrasound. Master’s thesis, Univ. B.C., Vancouver, 2004. [136] E. Turgay, S. E. Salcudean, and R. N. Rohling. Identifying mechanical properties of tissue by ultrasound. Ultrasound Med. Biol., 32:221–235, 2006. [137] I. Tutar, L. Gong, and S. Narayanan. Seed-based ultrasound and fluoroscopy registration using iterative optimal assignment for intraoperative prostate brachytherapy dosimetry. In Proc. SPIE, volume 6509, pages 14–21, 2007. [138] I. B. Tutar, R. Managuli, V. Shamdasani, P. S. Cho, S. D. Pathak, and Y. Kim. Tomosynthesis-based localization of radioactive seeds in prostate brachytherapy. Med. Phys., 30:3135–3142, 2003. [139] T. Varghese, E. E. Konofagou, J. Ophir, and et al. Direct strain estimation in elastography using spectral cross-correlation. Ultrasound Med. Biol., 26:1525–1537, 2000. [140] T. Varghese and J. Ophir. Performance optimization in elastography: Multi-compression with temporal stretching. Ultrason. Imaging, 18:193–214, 1996. [141] T. Varghese and J. Ophir. A theoretical framework for performance characterization of elastography: The strain filter. IEEE Trans. Ultrason., Ferroelect., Freq. Contr., 44:164– 172, 1997. [142] T. Varghese and J. Ophir. Analysis of elastographic contrast-to-noise ratio. Ultrasound Med. Biol., 24(6):915–924, 1998. [143] F. Viola and W. F. Walker. Imaging viscoelastic properties of the vitreous. In Proc. IEEE Ultrasonics Symp., pages 1623–1626, 2001. [144] F. Viola and W. F. Walker. A comparison of the performance of time-delay estimators in medical ultrasound. IEEE Trans. Ultrason., Ferroelect., Freq. Contr., 50:392–401, 2003. [145] Z. Wei, L. Gardi, D. B. Downey, and A. Fenster. Automated localization of implanted seeds in 3D TRUS images used for prostate brachytherapy. Med. Phys., 33:2404–2417, 2006. [146] Z. Wei, L. Gardi, D. B. Downey, and A. Fenster. Segmentation of implanted radioactive seeds in 3D TRUS images for intraoperative prostate brachytherapy. In Proc. IEEE ISBI, pages 876–879, 2006. [147] X. Wen and S. E. Salcudean. Enhancement of bone surface visualization using ultrasound radio frequency signals. In Proc. IEEE Ultrasonics Symp., 2007. 111  [148] J. Xue, F. Waterman, J. Handler, and E. Gressen. Localization of linked 125 I seeds in postimplant TRUS images for prostate brachytherapy dosimetry. Int. J. Radiat. Oncol. Biol. Phys., 62:912–919, 2005. [149] Y. Yu, S. T. Acton, and K. Thornton. Detection of radioactive seeds in ultrasound images of the prostate. In Proc. IEEE Int. Conf. Image Process., volume 2, pages 319–322, 2001. [150] R. Zahiri-Azar. Real-time imaging of elastic properties of soft tissue with ultrasound. Master’s thesis, Univ. B.C., Vancouver, 2005. [151] R. Zahiri-Azar and S. E. Salcudean. Motion estimation in ultrasound images using time domain cross correlation with prior estimates. IEEE Trans. Biomed. Eng., 53:1990–2000, 2006. [152] M. Zelefsky, Y. Yamada, C. Marion, and et al. Improved conformality and decreased toxicity with intraoperative computer-optimized transperineal ultrasound-guided prostate brachytherapy. Int. J. Radiat. Oncol. Biol. Phys., 55:956–963, 2003. [153] Y. Zhan and D. Shen. Automated segmentation of 3-D ultrasound prostate images using statistical texture-based matching method. In Proc. MICCAI, pages 688–696, 2003. [154] M. Zhang, M. Zaider, M. Worman, and G. Cohen. On the question of 3D seed reconstruction in prostate brachytherapy: The determination of X-ray source and film locations. Phys. Med. Biol., 49:335–345, 2004. [155] A. L. Zietman. Localized prostate cancer: Brachytherapy. Curr. Treat. Options Oncol., 3:429–436, 2002.  112  Appendix A  Ethics Certificate This Appendix presents the photocopy of the certificate of full board approval by the Clinical Research Ethics Board of the University of British Columbia.  113  114  Appendix B  Consent Form for Vibro-Elastography Volume Study This Appendix presents the photocopy of the subject information and consent form for the vibroelastography imaging of the prostate region of patients during the volume study, in the British Columbia Cancer Agency, Vancouver Centre.  115  116  117  118  119  Appendix C  Consent Form for Vibro-Elastography Study in the OR This Appendix presents the photocopy of the subject information and consent form for the vibroelastography imaging of the prostate region of patients in the operating room pre-operatively, in the British Columbia Cancer Agency, Vancouver Centre, and for the MRI scan of the prostate region in the UBC hospital.  120  121  122  123  124  125  126  Appendix D  Data Acquisition Protocols D.1  Data Acquisition Protocol for Vibro-Elastography Volume Study  D.1.1  Purpose  We wish to assess the visibility of the prostate with real-time ultrasound elastography. This is a pilot study with a limited number of patients (at most 20) with the goal of assessing and tuning the imaging equipment before carrying out a more extensive study involving intra-operative data acquisition.  D.1.2  Study location and duration  The data acquisition will take place in the same room where standard pre-operative volume studies are carried out, immediately after the standard clinical work will be completed. The additional examination time will be limited to 10 minutes.  D.1.3  Summary of the procedure and equipment  An additional elastographic prostate volume study will be carried out after the standard preoperative ultrasound volume study. For this study, a different ultrasound machine, TRUS probe and brachytherapy stepper will be used. The prostate localization will be carried out in the usual manner. After the prostate is localized, a sagittal scan will be taken. During this scan, a small vibration will be applied to the ultrasound transducer by a small motor. Ultrasound strain images will be acquired. The following equipment will be moved to the pre-op scan room prior to the patient’s arrival: 1. The Ultrasonix ultrasound machine with a dual-plane endorectal TRUS probe. The Ultrasonix will be placed adjacent to the Siemens machine. 2. The CIVCO EXII stepper with the TRUS cradle and motorized shaker mounted on it. This stepper will replace the CIVCO Classic stepper you normally use for volume scans. 3. The TRUS depth and rotation angles monitor. This is an electronics box the size of a large Kleenex tissue box that has a small display showing the TRUS insertion depth in millimeters and the TRUS angle in degrees. It will be placed on the Ultrasonix machine where it can be seen and accessed easily. 4. Power supplies and cables for the motorized shaker and the TRUS depth and angle monitor.  127  D.1.4  Suggested equipment setup  1. Place the above-mentioned equipment in the pre-op scan room (Xu and Staff) 2. Connect the TRUS monitor to the Ultrasonix machine (Xu) 3. Connect the shaker with the Ultrasonix machine (Xu) 4. Start running the Ultrasonix machine software and ensure its operation (Xu) 5. Start running the shaker software and ensure its operation (Xu) 6. Put a condom onto the TRUS probe with gels (Staff) 7. Place the TRUS probe on the EXII Brachystand (Radiation oncologist performing volume study)  D.1.5  Data Acquisition Steps  1. Probe switch: Replace the Siemens probe and CIVCO Classic stepper with the Ultrasonix probe and the CIVCO EXII stepper (Staff or Xu). 2. Probe insertion: Insert the probe to the patient’s rectum, adjust the position of the TRUS probe, and locate the base plane of the prostate in the usual manner (radiation oncologist). Record the base and apex positions and reset the encoder at the base (Xu). Note that the Classic stepper and the EXII stepper are mounted on the same CIVCO Microtouch stabilizer with the same type of mount. Therefore the transducer adjustment required should be minimal (we hope for less than 2 minutes). 3. Vibro-elastography tuning: Move the probe to the base plane (radiation oncologist). Activate the sagittal scan and adjust the amplitude and frequency for good strain image contrast (radiation oncologist and Xu) (2 minutes). 4. Vibro-elastography data acquisition. Rotate the cradle (probe) clockwise to -45 deg (radiation oncologist). From -45 deg, rotate the cradle (radiation oncologist) counterclockwise and stop every 5 deg (radiation oncologist) when prompted by Xu, to +45 deg. Between -5 deg and 5 deg, stop every 1 deg (radiation oncologist). At each angle, acquire data for 10 - 15 seconds (Xu). (4-5 minutes). 5. Additional imaging. Optional. Would like to see if elastography can provide improved imaging of the nerves and vasculature adjacent to the prostate (angles required may be larger than 45 deg). (radiation oncologist with Xu’s help) (2 minutes). 6. Probe retraction: Retract the probe from the patient’s rectum (radiation oncologist), and remove the condom (Staff). Disconnect the probe from the Ultrasonix machine and soak it (Xu). Disconnect the shaker and encoder monitor from the Ultrasonix machine and the Microtouch brachystand (Xu).  128  D.2  Data Acquisition Protocol for Vibro-Elastography Study in the OR  D.2.1  Study location and duration  The data acquisition will take place in the OR where standard prostate brachytherapy is carried out, before the standard clinical work starts and after the standard clinical work is completed. The additional examination time will be limited to 15 minutes.  D.2.2  Equipment  The following equipment will be moved to the OR prior to the patient’s arrival: 1. The Ultrasonix ultrasound machine with a dual-plane endorectal TRUS probe. The Ultrasonix will be placed adjacent to the B-K machine. 2. The CIVCO EXII stepper with the TRUS cradle and motorized shaker mounted on it. This stepper will replace the CIVCO Classic stepper you normally use. 3. Power supplies and cables for the motorized shaker. 4. Digital protractor.  D.2.3  Suggested equipment setup  1. Place the above-mentioned equipment in the OR (Shawn) 2. Connect the shaker to the Ultrasonix machine (Shawn) 3. Start running the Ultrasonix machine software and ensure its operation (Shawn) 4. Start running the shaker software, calibrate the shaker (“Calibrate” - “Execute roll” “Home”), and initialize the band-pass vibration (Shawn). 4. Put a condom onto the TRUS probe with gels (Staff) 5. Place the TRUS probe on the EXII Brachystand (Shawn)  D.2.4  Data Acquisition Steps  Before the brachytherapy procedure 1. Adjust the position and height of the patient table and C-arm to make sure that there is sufficient space for the C-arm to rotate between -10 and 10 deg (radiation oncologist, staff, and Shawn) without collision. Measure the flip-flop angle of the C-arm (Shawn). (3 minutes) 2. Probe insertion: Insert the probe to the patient’s rectum, adjust the position of the TRUS probe, and locate the base plane of the prostate in the usual manner (radiation oncologist). 3. Vibro-elastography tuning: Move the probe to the base plane (radiation oncologist). Activate the sagittal scan, adjust the probe position and imaging depth to make sure that we can see the whole prostate, and check the image quality at 0 and/or 45 deg (radiation oncologist and Shawn). 4. Volume study - elastography: Start running the software (roll with band-pass vibration and increments) and collect ultrasound data (Shawn). (2 minutes) 129  5. Volume study - no vibration (optional): Start running the software (roll without vibration and increments) and collect ultrasound data (Shawn). (1 minute) 6. Probe switch: Replace the Ultrasonix probe and the CIVCO EXII stepper with the B-K probe and CIVCO Classic stepper (staff and Shawn). After the brachytherapy procedure 7. Probe switch: Replace the B-K probe and CIVCO Classic stepper with the Ultrasonix probe and the CIVCO EXII stepper (staff and Shawn). Measure the angle of the stepper (Xu). 8. Probe insertion: Insert the probe to the patient’s rectum, adjust the position of the TRUS probe, and locate the base plane of the prostate in the usual manner (radiation oncologist). 9. Volume study - seeds: Move the probe to the base plane (radiation oncologist). Activate the sagittal scan, and adjust the probe position to make sure that we see can the whole prostate (radiation oncologist). Start running the software (roll without vibration and increments) and collect ultrasound data (Xu). (1 minute) 10. Volume study - vibration (optional): Click the ”Home” button to return the probe back to 0 deg (Xu). Start running the software (roll with single-frequency vibration and increments) and collect ultrasound data (Xu). (2 minute) 11. Rotate the cradle (probe) back to 0 deg (radiation oncologist). Measure the distance from the X-ray source to the probe (Xu). 12. Fluoro scan: Retract the probe, expose and save the fluoro image without the probe (staff or Xu). Then insert the probe to the rectum to make sure that there is no overlap between the probe and the seeds (radiation oncologist, staff and Shawn). Expose and save the fluoro image (staff or Shawn). Rotate the C-arm to -10 deg, expose and save the fluoro image (staff or Xu). Then rotate the C-arm to -5, 0, +5, and +10 deg, respectively, expose and save the fluoro images (staff or Xu), and record the readings of the digital protractor for all the fluoro images (Xu and Nick). (3 minutes)  130  Appendix E  Communication Protocol between the Data Acquisition Program and VE Actuation System The real-time data acquisition program and the VE actuation system program are running on the same computer (the Sonix RP machine in the experiments) independently. The communication between the two programs has been implemented by two shared text files: one for rotation angles of the cradle and transducer, and the other for vibrating positions of the cradle and transducer. The two files are read and written independently as well. The paths of the two shared files are configurable in the config files associated with the two programs, respectively.  E.1  Writing Shared Files  When the cradle and transducer are automatically rotated from the start angle to the stop angle, the most recently sampled value of the rotation angle is overwritten to the rotation angle file, with resolution of 0.1◦ . When the cradle and probe are rotated from the zero-position (0◦ ) to the start angle, stopped at the stop angle or stopped anytime in the middle, or rotated from the last angle back to the zero-position, the angle value overwritten to the file is -9999, instead of the actual rotation angles. The maximal angle reading range is from -48.0 to 50.0. The file is locked when it is written. Due to the electrical performance limitation, the maximal sampling rate of rotation angles is 50 Hz currently. When the cradle and transducer are vibrated, the most recently sampled value of the vibrating position is overwritten to the vibrating position file. The sampled value is provided by the encode with the vibrating motor. The relationship between the sampled value and the actual transducer position has not been accurately determined yet. The file is locked when it is written. Due to the electrical performance limitation, the maximal sampling rate of vibrating positions is 50 Hz currently.  E.2  Reading Shared Files  While the real-time data acquisition program is running, it keeps reading the two shared text files every time a new frame of RF data is acquired. The reading rate is equal to the data frame rate and is typically lower than the writing rate. When the rotation angle reading is changed from -9999 to any other values, RF or B-mode 131  data acquisition will be automatically started. The RF and/or B-mode data are stored in memory. The rotation angle and shaking position of the transducer are saved in the head of each RF or B-mode data frame, in the data type of float. When the angle reading is changed from any other values to -9999, data collection will be automatically stopped if it has been started. Then the RF and/or B-mode data in memory are dumped to the hard drive.  132  Appendix F  Recipes for Making Phantoms F.1  Equipment  The following equipment is used for making phantoms: a beaker, a hot stage, a scale or measuring spoon, a stirrer, and a thermometer. Additionally, agar or gelatin phantoms should be kept in a container in a fridge and covered by food wrappers if possible.  F.2  Making agar phantoms  An agar phantom is typically made in the following steps: 1. Weigh the agar and/or cellulose powders by a scale or measuring spoon, according to your requirements. Typically the percent of agar in mass is 2% - 3%, and the percent of agar in mass is 1% - 3%. 2. Pour water (or deionized water) in a beaker to the expected volume, according to the requirements 3. Heat the water (in the beaker) on a hot stage to 90 - 100 ◦ C (lower than 100 ◦ C) 4. Put the agar powders into the water, and stir the mixture at the same time 5. Heat and stir the mixture for a few minutes on the hot stage, and then put the cellulose powders into the water as scatterers if needed 6. Stop heating, wait for the mixture to cool down, and keep stirring the water from time to time 7. Pour the mixture into a container when the temperature drops below 60 ◦ C, and then leave the container in a fridge  F.3  Making gelatin phantoms  An gelatin phantom is typically made in the following steps: 1. Weigh the gelatin and/or cellulose powders by a scale or measuring spoon, according to your requirements. Typically the percent of gelatin in mass is 10% - 15%, and the percent of agar in mass is 1% - 3%. 2. Pour water (or deionized water) in a beaker to the expected volume, according to the requirements 3. Heat the water (in the beaker) on a hot stage to 40 - 50 ◦ C 4. Put the gelatin powders into the water, and stir the mixture at the same time 133  5. Heat and stir the mixture for a few minutes on the hot stage, and then put the cellulose powders into the water as scatterers if needed 6. Stop heating, wait for the mixture to cool down, and keep stirring the water from time to time 7. Pour the mixture into a container when the temperature drops between 28 - 29 ◦ C (not lower than 28 ◦ C), and then leave the container in a fridge  134  

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.24.1-0071078/manifest

Comment

Related Items