Information Fusion for ProstateBrachytherapy PlanningbySaman NouranianB.Sc., Amirkabir University of Technology, 1999M.Sc., University of Tehran, 2003A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Electrical and Computer Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)June 2016c© Saman Nouranian 2016AbstractLow-dose-rate prostate brachytherapy is a minimally invasive treatment ap-proach for localized prostate cancer. It takes place in one session by perma-nent implantation of several small radio-active seeds inside and adjacent tothe prostate. The current procedure at the majority of institutions requiresplanning of seed locations prior to implantation from transrectal ultrasound(TRUS) images acquired weeks in advance. The planning is based on a setof contours representing the clinical target volume (CTV). Seeds are man-ually placed with respect to a planning target volume (PTV), which is ananisotropic dilation of the CTV, followed by dosimetry analysis. The mainobjective of the plan is to meet clinical guidelines in terms of recommendeddosimetry by covering the entire PTV with the placement of seeds. The cur-rent planning process is manual, hence highly subjective, and can potentiallycontribute to the rate and type of treatment related morbidity.The goal of this thesis is to reduce subjectivity in prostate brachytherapyplanning. To this end, we developed and evaluated several frameworks toautomate various components of the current prostate brachytherapy plan-ning process. This involved development of techniques with which targetvolume labels can be automatically delineated from TRUS images. A seedarrangement planning approach was developed by distributing seeds withrespect to priors and optimizing the arrangement according to the clinicalguidelines. The design of the proposed frameworks involved the introductionand assessment of data fusion techniques that aim to extract joint informa-tion in retrospective clinical plans, containing the TRUS volume, the CTV,the PTV and the seed arrangement. We evaluated the proposed techniquesusing data obtained in a cohort of 590 brachytherapy treatment cases fromthe Vancouver Cancer Centre, and compare the automation results with theclinical gold-standards and previously delivered plans. Our results demon-strate that data fusion techniques have the potential to enable automaticplanning of prostate brachytherapy.iiPrefaceThis thesis is primarily based on two journal publications and five confer-ence papers, resulting from the collaboration between multiple researchers.The author was responsible for development, implementation and evalua-tion of the method and the production of the manuscripts. All co-authorshave contributed to the editing of the manuscripts and providing feedbackand comments. The research conducted in this study was undertaken underthe approval of the UBC BCCA Research Ethics Board, certificate numberH13-01983 titled ”Information Fusion Technology to Optimize Brachyther-apy Procedure”.The study from Chapter 2 is presented at:• S. Nouranian, S. S. Mahdavi, I. Spadinger, W. J. Morris, S. E. Salcud-ean, and P. Abolmaesumi, Multi-atlas-based Automatic 3D Segmenta-tion for Prostate Brachytherapy in Transrectal Ultrasound Images, inProceedings of SPIE Medical Imaging, 2013, vol. 8671, pp. 86710O1–7.• S. Nouranian, S. S. Mahdavi, I. Spadinger, W. J. Morris, S. E. Sal-cudean, and P. Abolmaesumi, An Automatic Multi-atlas Segmenta-tion of the Prostate in Transrectal Ultrasound Images Using PairwiseAtlas Shape Similarity, in Medical Image Computing and Computer-Assisted Intervention: MICCAI, 2013, vol. 8150, pp. 173-180.An extended version of the proposed methodology from Chapter 2 validatedon a larger dataset has been published in:• S. Nouranian, S. S. Mahdavi, I. Spadinger, W. Morris, S. Salcudean,and P. Abolmaesumi, A Multi-Atlas-Based Segmentation Frameworkfor Prostate Brachytherapy, IEEE Trans. Med. Imaging, vol. 34, no.4, pp. 950-961, 2015.Dr. Mahdavi provided a dataset of previously studied TRUS volumes andsegmentation results of the semi-automatic approach part of the currentclinical software package. Scientific inputs of Prof. Abolmaesumi and Prof.iiiPrefaceSalcudean helped with development and implementation of the methodol-ogy. Dr. Spadinger and Dr. Morris provided advice on clinical evaluationapproaches. All authors helped to improve the manuscript structure by theirvaluable comments and feedback.Studies presented in Chapter 3 are presented at:• S. Nouranian, M. Ramezani, S. S. Mahdavi, I. Spadinger, W. J. Morris,S. E. Salcudean, and P. Abolmaesumi, Fast Prostate Segmentationfor Brachytherapy based on Joint Fusion of Images and Labels, inProceedings of SPIE Medical Imaging, 2014, vol. 9036, pp. 90361A1-7.• S. Nouranian, M. Ramezani, S. S. Mahdavi, I. Spadinger, W. J. Mor-ris, S. E. Salcudean, and P. Abolmaesumi, Data Fusion for PlanningTarget Volume and Isodose Prediction in Prostate Brachytherapy, inProceedings of SPIE Medical Imaging, vol. 9415. pp. 94151I1-7, 2015.Dr. Mahdavi provided the dataset used in evaluation of the presentedmethod in these two papers. Prof. Abolmaesumi, Prof. Salcudean andDr. Ramezani have contributed to the development and improvement ofthe methodology. Dr. Spadinger and Dr. Morris provided advice on clin-ical evaluation approaches. All authors helped to improve the manuscriptstructure by their valuable comments and feedback.A study described in Chapter 4 has been published in:• S. Nouranian, M. Ramezani, I. Spadinger, W. Morris, S. Salcudean,and P. Abolmaesumi, Learning-based Multi-label Segmentation of Tran-srectal Ultrasound Images for Prostate Brachytherapy, IEEE Trans.Med. Imaging, vol. 35, no. 3, pp. 921-932, 2016.Prof. Abolmaesumi, Prof. Salcudean and Dr. Ramezani provided scientificinput for development and improvement of the methodology. Dr. Spadingerand Dr. Morris provided advice on clinical evaluation approaches. All au-thors helped to improve the manuscript structure by their valuable com-ments and feedback.A study described in Chapter 5 is presented at:• S. Nouranian, M. Ramezani, I. Spadinger, W. J. Morris, S. E. Salcud-ean, and P. Abolmaesumi, Automatic Prostate Brachytherapy Pre-planning Using Joint Sparse Analysis, in Medical Image Computingand Computer-Assisted Intervention: MICCAI, 2015, vol. 9350, pp.415–423.ivPrefaceProf. Abolmaesumi, Prof. Salcudean and Dr. Ramezani provided scientificinput for development and improvement of the methodology. Dr. Spadingerand Dr. Morris provided advice on clinical evaluation approaches. All au-thors helped to improve the manuscript structure by their valuable com-ments and feedback.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xivDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Brachytherapy . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Uncertainties in LDR Brachytherapy . . . . . . . . . . . . . . 31.3 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.3.1 Target Volume Segmentation . . . . . . . . . . . . . . 81.3.2 Seed Arrangement Planning . . . . . . . . . . . . . . . 101.4 Proposed Solution . . . . . . . . . . . . . . . . . . . . . . . . 101.4.1 Objective . . . . . . . . . . . . . . . . . . . . . . . . . 111.4.2 Contributions . . . . . . . . . . . . . . . . . . . . . . . 121.5 Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121.6 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . 142 Atlas-based Fusion for Brachytherapy Preplanning . . . . . 172.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.2.1 Target-specific Atlas Selection . . . . . . . . . . . . . . 192.2.2 Filtered Multi-atlas Fusion . . . . . . . . . . . . . . . 21viTable of Contents2.2.3 Clinical Target Volume Estimation . . . . . . . . . . . 222.3 Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232.3.1 Segmentation Evaluation . . . . . . . . . . . . . . . . 252.4 Experiments and Results . . . . . . . . . . . . . . . . . . . . . 272.4.1 Mean Atlas Generation . . . . . . . . . . . . . . . . . 282.4.2 Global and Sector-based Performance Measures . . . . 282.4.3 Computational Cost . . . . . . . . . . . . . . . . . . . 322.4.4 Comparison of Different Fusion Methods . . . . . . . . 332.5 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . 352.5.1 Comparison with Other Methods . . . . . . . . . . . . 362.5.2 Dataset Dependency . . . . . . . . . . . . . . . . . . . 373 ICA-based Fusion for Brachytherapy Preplanning . . . . . 393.1 Clinical Target Volume Estimation from TRUS Images . . . . 393.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . 393.1.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 403.1.3 Results and Discussion . . . . . . . . . . . . . . . . . . 433.1.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 453.2 Planning Target Volume and Minimum Prescribed IsodoseEstimation from TRUS Images . . . . . . . . . . . . . . . . . 463.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . 463.2.2 Method . . . . . . . . . . . . . . . . . . . . . . . . . . 473.2.3 Experiments . . . . . . . . . . . . . . . . . . . . . . . 513.2.4 Results and Discussion . . . . . . . . . . . . . . . . . . 513.2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 534 Sparsity-based Fusion for Brachytherapy Preplanning . . . 554.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 584.2.1 Problem Definition and Observation Space Formation 584.2.2 Joint Sparse Dictionary Learning . . . . . . . . . . . . 604.2.3 Multi-label Estimation . . . . . . . . . . . . . . . . . . 614.3 Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.4 Experiments and Results . . . . . . . . . . . . . . . . . . . . . 634.4.1 Model Parameter Tuning . . . . . . . . . . . . . . . . 634.4.2 Multi-label Estimation . . . . . . . . . . . . . . . . . . 644.4.3 Dataset Adaptability Analysis . . . . . . . . . . . . . 664.4.4 Comparison of the CTV Labeling with Other Methods 674.4.5 Multi-label Model Comparison with Single-label Model 684.4.6 Execution Time . . . . . . . . . . . . . . . . . . . . . . 69viiTable of Contents4.5 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . 704.5.1 Clinical Significance . . . . . . . . . . . . . . . . . . . 704.5.2 Computational Cost . . . . . . . . . . . . . . . . . . . 714.5.3 Dataset Adaptability . . . . . . . . . . . . . . . . . . . 714.5.4 Limitations . . . . . . . . . . . . . . . . . . . . . . . . 725 Automatic Seed Plan Estimation . . . . . . . . . . . . . . . . 785.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 785.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 795.2.1 Joint Sparse Dictionary Learning . . . . . . . . . . . . 795.2.2 Seed Plan Estimation . . . . . . . . . . . . . . . . . . 805.2.3 Plan Optimization . . . . . . . . . . . . . . . . . . . . 815.3 Experiments and Results . . . . . . . . . . . . . . . . . . . . . 875.4 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . 886 Conclusion and Future Work . . . . . . . . . . . . . . . . . . . 916.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 93Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98viiiList of Tables2.1 Volumetric, contour and surface error statistics (mean±standarddeviation) calculated from 10 random dataset partitioningsare shown for the proposed automatic segmentation approachand compared with the approach of Mahdavi et al. [1]. . . . . 292.2 Grand mean and pooled standard deviation of Dice (DSC)calculated for 10 different trials of randomly dividing datainto 60-20-20% training-validation-test sets. . . . . . . . . . . 322.3 Grand mean and pooled standard deviation of mean surfacedistance (MSD) calculated for 10 different trials of randomlydividing data into 60-20-20% training-validation-test sets. . . 333.1 Comparison of various evaluation metrics between the pro-posed segmentation approach and the manually segmentedcontours determined by an expert clinician. . . . . . . . . . . 443.2 Grand mean and pooled standard deviation of Verr, Vdiff andV 100 calculated for the estimated PTV and mPD volume in10 different random training/testing datasets. . . . . . . . . . 534.1 Mean and standard deviation of Verr calculated for nine anatom-ical sectors of the prostate over CTV and PTV. . . . . . . . . 654.2 Comparison to two other brachytherapy CTV segmentationmethods for sparse TRUS volumes. “?” and “•” indicate astatistically significant difference between mean and standarddeviation of the error measured by the proposed multi-labelapproach and each of the reported methods using paired t-test and Levene’s test, respectively. . . . . . . . . . . . . . . . 675.1 Clinical guidelines at our local cancer centre: . . . . . . . . . 825.2 Mean and standard deviation of the plan quality metrics cal-culated for actual, estimated and random plans with andwithout optimization. . . . . . . . . . . . . . . . . . . . . . . 88ixList of Figures1.1 Prostate cancer treatment methods. . . . . . . . . . . . . . . 21.2 Volume study in brachytherapy . . . . . . . . . . . . . . . . . 31.3 Clinical and planning target volumes . . . . . . . . . . . . . . 41.4 Schematic of a prostate brachytherapy implantation procedure 51.5 Current prostate brachytherapy preplanning requires user in-tervention as shown by circles. . . . . . . . . . . . . . . . . . 111.6 Recovering ultrasound images from retrospective dataset ofpatient charts . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.1 Block diagram of the proposed MAS-based TRUS segmenta-tion framework. . . . . . . . . . . . . . . . . . . . . . . . . . . 202.2 Process of making TRUS contours smooth and symmetric . . 242.3 Prostate sectioning for segmentation evaluation . . . . . . . . 262.4 Correlation analysis between registration performance andsegmentation overlap . . . . . . . . . . . . . . . . . . . . . . . 282.5 Atlas at the mid-gland slice . . . . . . . . . . . . . . . . . . . 292.6 DSC values computed for nine anatomical sectors of the prostate. 302.7 Average operating computation time with total number of100 atlases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 302.8 Comparison between different fusion algorithms . . . . . . . . 312.9 Two examples of the automatic multi-atlas based segmenta-tion result. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 342.10 Performance comparison between the proposed automatic ap-proach and the semi-automatic method . . . . . . . . . . . . 373.1 Block diagram of the jICA based segmentation algorithm . . 423.2 Contour smoothing process . . . . . . . . . . . . . . . . . . . 423.3 Contours for six equally-spaced slices (base to apex) of twotypical segmentation results. . . . . . . . . . . . . . . . . . . . 443.4 Histogram of the volumetric error and surface distance measure 453.5 Isodose prediction pipeline . . . . . . . . . . . . . . . . . . . . 473.6 Reconstruction error calculated for PTV and Isodose . . . . . 52xList of Figures3.7 Verr statistics shown with boxplots from 10 random datasetpartitioning into train and test sets. . . . . . . . . . . . . . . 523.8 Gold standard and predicted TRUS, PTV and mPD volume . 544.1 Sample CTV and PTV contours used for treatment of twopatients. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 564.2 Illustration of the proposed system for simultaneous auto-matic labeling of the TRUS volumes. . . . . . . . . . . . . . . 594.3 Hyper-parameter tuning by grid search on sparsity level anddictionary size . . . . . . . . . . . . . . . . . . . . . . . . . . 644.4 Visualized sparse dictionary atoms . . . . . . . . . . . . . . . 654.5 Label estimation errors in different anatomical sectors of theprostate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 664.6 Two typical results of multi-label estimation from TRUS vol-ume. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.7 Learning-based multi-label segmentation adaptability analysis. 754.8 Cumulative histogram of the CTV estimation . . . . . . . . . 764.9 Comparison of single-label models with the multi-label approach 775.1 Block diagram of the proposed seed planning framework . . . 795.2 Initial plan estimation using the joint sparse dictionary learn-ing algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 815.3 Architecture of proposed optimization algorithm based onsimulated annealing. . . . . . . . . . . . . . . . . . . . . . . . 835.4 Seed arrangement correction and annealing process for oneiteration of simulated annealing . . . . . . . . . . . . . . . . . 855.5 Seed rearrangement correction based on calculated probabil-ity values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 855.6 Comparison between statistics of plan quality metrics . . . . 89xiGlossaryAAM Active Appearance Model.ABS American Brachytherapy Society.BCCA British Columbia Cancer Agency.CT Computed Tomography.CTV Clinical Target Volume.DRE Digital Rectal Examination.DSC Dice Similarity Coefficient.GPU Graphical Processing Unit.HD Hausdorff Distance.ICA Independent Component Analysis.ICRU International Commission of Radiation Units.jICA joint Independent Component Analysis.K-SVD K-Singular Value Decomposition.KSVM Kernel Support Vector Machines.LDR Low Dose Rate.LOP Logarithmic Opinion Pool.LWMV Locally Weighted Majority Voting.xiiGlossaryMAD Mean Absolute Radial Distance.MAP Maximum A Posteriori.MAS Multi-Atlas Segmentation.MAXD Maximum Absolute Radial Distance.mPD Minimum Prescribed Dose.MR Magnetic Resonance.MSD Mean Surface Distance.OR Operating Room.PCA Principal Component Analysis.PSA Prostate Specific Antigen.PTV Planning Target Volume.RBM Restricted Bultzman Machines.SA Simulated Annealing.SIMPLE Selective and Iterative Method for Performance Level Estima-tion.SSD Sum of Squared Differences.SSM Statistical Shape Model.STAPLE Simultaneous Truth and Performance Level Estimation.TRUS Transrectal Ultrasound.VCC Vancouver Cancer Centre.xiiiAcknowledgementsFirst and foremost, I would like to express my sincerest gratitude to Pro-fessor Purang Abolmaesumi, my advisor, for the opportunities you gave tome. Thank you for your trust, patience and invaluable support that openedthis chapter of my life. Thanks for the generosity, advice and friendship youoffered me throughout my experience at UBC.I am very grateful for being associated with Drs. James Morris andIngrid Spadinger, whose insightful advices helped me throughout the courseof this endeavor. Their patience, dedication and sincere advice in everysingle step of the thesis have been truly crucial. Thanks for being alwayssupportive and enthusiastic.I would like to thank supervisory committee Professors Salcudean andRohling for their valuable comments and advices.I would like to particularly thank Dr. Mahdi Ramezani for his not onlyscientific support but also friendship and enthusiasm shared with me duringyears I started at UBC.I would like to thank my dearest friend Dr. Mehdi Moradi for his un-conditional support and friendship since I started my academic journey.My sincerest thanks to all current and past RCL members at UBC, Vickie,Abtin, Farhad, Siavash, Alex, Jeff, Parmida, Shekoofeh, Rohit, Greg, Philipand many others.I owe immensely my true friend and wife Samira for her invaluable sup-port in all small and big steps I’ve taken during this endeavor. And deepmy heart, I thank my mom, dad and brother who constantly supported andinspired me to start this journey. You will never be forgotten.xivDedicationI dedicate this to my father who left fingerprints of grace on my life; mymother and my wife without whom none of my success would be possible ...xvChapter 1Introduction1.1 BrachytherapyProstate cancer is one of the frequently diagnosed types of cancer in NorthAmerica. According to American and Canadian Cancer Societies, it hasaffected the life of at least one in six men in North America in 2013 [2, 3].Recent advances in diagnosis and treatment methods proposed for prostatecancer have shown a significant increase in the rate of survival; e.g. theCanadian Cancer Society reported five-year relative survival ratio of morethan 95% for patients diagnosed between 2006 and 2008 comparing to the86% for the range of 1992 to 1994 [3]. Prostate cancer is routinely detectedby Digital Rectal Examination (DRE) and measurement of the ProstateSpecific Antigen (PSA). Although PSA level is often elevated in the presenceof cancer, reliable diagnosis leading to staging of the cancer is usually madeby pathological analysis of the sample tissues collected from the prostatecore using biopsy needles under guidance of transrectal ultrasound (TRUS)images.Several treatment methods have been proposed to eliminate or removethe cancerous tissue or the entire gland, while minimizing unnecessary mor-bidity and better preserving patients’ quality of life (see Figure 1.1). Low-dose-rate (LDR) prostate brachytherapy is an effective treatment option forlocalized prostate cancer with low and intermediate risk disease (based onNational Comprehensive Cancer Network1 guidelines). It is associated witha good long-term survival rate [4, 5], requires less patient hospitalizationwhen compared to radical prostatectomy, and is delivered in one session asopposed to alternative methods such as external beam therapy that need tobe administered over a period of multiple sessions.In LDR brachytherapy, small radioactive seeds are loaded as multi-seed”trains” into several needles in predefined patterns. They are permanentlyimplanted through the perineum into the prostate and its adjacent tissueusing a grid template (see Figure 1.4). Throughout this thesis, we refer to1http://www.nccn.org11.1. Brachytherapythe LDR brachytherapy treatment approach that has been developed in theBritish Columbia Cancer Agency (BCCA) [6]. This program is founded in1997 by radiation oncologists and medical physicists at the BCCA. Treat-ment protocol and guidelines are a combination of an in-house manual seedplanning algorithm and evolved version of the Seattle preplanning expe-rience [7]. According to the BCCA guidelines, a typical brachytherapyprocedure can be divided into four main steps, usually taking place daysapart; volume study, seed planning, implantation, and post-implant analysis.Prostate Cancer Treatment Active Surveillance Whole-gland Ablation LDR Brachytherapy HDR Brachytherapy External Beam Photons Proton beam radiation therapy Radical Prostatectomy Cryosurgery Systemic Therapies Hormone therapy Biologic therapy Chemotherapy Focal Therapy HIFU Focal cryoablation Figure 1.1: Prostate cancer treatment methods.The volume study step involves determination of the Clinical Target Vol-ume (CTV) or the target eradicable anatomy boundary. Standard-of-carein this process is to use an ultrasound machine and a transrectal ultrasoundprobe to acquire a sparse volumetric representation of the gland by a set of5 mm apart 2-D parallel images from the base (bladder) to the apex (pelvicfloor). Subsequently, all these 2-D images are contoured by an expert ra-diation oncologist to provide the CTV, which closely but not necessarilyexactly follows the prostate boundary (Figure 1.2). CTV cannot be accu-rately defined for an individual patient; but is meant to contain all cancerousanatomy that requires to be adequately treated.In seed planning step, the number and distribution of seeds and nee-dles is determined based on some dosimetry analysis. First, margins areadded to the CTV, to obtain the Planning Target Volume (PTV). PTV isdetermined to account for uncertainties of delivering the planned dose tothe CTV. Generally, the dilation parameters are suggested in clinical guide-21.2. Uncertainties in LDR Brachytherapy Base Apex TRUS probe Apex Base TRUS probe Figure 1.2: Volume study in brachytherapy requires prostate boundary de-lineation in a series of equally-spaced parallel 2-D transrectal ultrasoundimages.lines. Figure 1.3 illustrates the relation between CTV and PTV based onthe guidelines used at the Vancouver Cancer Centre (VCC) and the BCCA.A predefined set of margins are used to spare the CTV volume by 0.5 cmsuperiorly, inferiorly and laterally, 0.3 cm anteriorly and 0.0 cm posteri-orly. After this anisotropic dilation of the CTV, the shape and size of thePTV may then be further altered slice by slice as needed to improve thecompromise between dose coverage and dose conformality. These guidelinesrecommend PTV determination, seed placement and needle configurationsto be symmetric across the patient mid-lobe.In the procedure day, patient is anaesthetised and positioned to the sameposture and condition as the volume study day. Radiation oncologist followsthe preplanned seed map to perform implantation and to deliver minimumprescribed dose (mPD) to the CTV accurately.1.2 Uncertainties in LDR BrachytherapyBrachytherapy procedure outcome and morbidity rate is assessed by pa-tient follow up records and recurrence state, months and years after treat-ment. Like other prostate cancer treatment methods, adjacency of urinary31.2. Uncertainties in LDR BrachytherapyApex+0.5 cm plane Axial View PTV seminal vesicles Base plane = Base+0.0 cm First plane with CTV = Base+0.5 cm Apex plane CTV Sagittal View Coronal View Figure 1.3: Planning target volume (PTV) is anisotropic dilation of clinicaltarget volume (CTV) to account for implantation variability.The significance of the intra-operative implantation uncertainty has ledthe American Brachytherapy Society (ABS) to recommend post-implantanalysis for every patient [8, 9]. CT imaging is unanimously recommendedas the standard protocol to detect seed positions in order to reconstruct thedose [10]. It is recommended that CT imaging should take place at approx-imately one month post-implant, at which time post-operative swelling haslargely resolved [11, 12]. Alternatively, CT imaging may be acquired imme-diately after the procedure, when there may be significant edema, in order toobtain more timely feedback on implant quality. As part of the standard-of-care at BCCA, day-0 CT is acquired right after the patient post-anaestheticrecovery.Implantation quality is routinely assessed on the basis of the post-implantD90 and V100 parameters. The D90 parameter is defined as the minimumdose in Gy received by 90% of the prostate volume and the V100 parameter is41.2. Uncertainties in LDR BrachytherapyFigure 1.4: Schematic of a prostate brachytherapy implantationprocedure2.percentage of the prostate volume receiving at least the minimum prescribeddose. Based on the criteria defined for these two parameters, implantationquality is classified as excellent, good or suboptimal.In addition to the post-implant analysis based on CT, brachytherapytreatment outcome is assessed by periodic monitoring of PSA and testos-terone levels over the course of several years after the procedure. Commonside effects such as rectal wall and urinary toxicity are also monitored.and erectile nerve bundles to the cancerous tissue causes prevalent post-procedural complications such as urinary incontinence, erectile dysfunctionand rectal wall toxicity. However, this rate of morbidity varies from patientto patient.Although some gene-related and physiological reasons play a key role indeveloping recurrent disease [13], brachytherapy outcome is affected by sev-eral obstacles to accurately deliver the intended dose to the target anatomy.Treatment-related sources of uncertainty can be divided into two groups:1) uncertainties associated with pre-operative tasks in determining the tar-get volume and the seed plan (preplanning); and 2) uncertainties stemmingfrom intra-operative tasks to deliver the plan.Prostate brachytherapy preplanning, as defined at the BCCA, involves2Anatomic schematic is retrieved from the book “Fast Facts: Prostate Cancer” (7thedition) by Kirby, Roger.51.2. Uncertainties in LDR BrachytherapyTRUS volume collection, CTV and PTV delineation, and determination ofthe seed arrangement w.r.t. the PTV. Accurate determination of the CTVis desirable in order to deliver sufficient radiation dose to the prostate whileminimizing dose to the urethra and the surrounding tissues, such as bladderand rectum. CTV segmentation error, in combination with other inherentsources of error (e.g. seed delivery error), can produce unnecessary morbid-ity, such as rectal wall toxicity, if the CTV is overestimated posteriorly, orlead to under-treatment, if the CTV is underestimated.CTV delineation is a cumbersome task that requires drawing of a con-tour for each 2-D TRUS image that mainly but not necessarily follows theprostate boundaries. Typically, a volume study image set consists of 7-14ultrasound images at 5 mm spaced axial planes. These images are oftenaffected by speckles, shadowing and reverberation artifacts, and the bound-ary of the prostate is not reliably visible, especially in the base and theapex. These characteristics of the TRUS images, in conjunction with thecurrent CTV delineation process which is either manual or semi-automatic,make delineation of the CTV a tedious task and vulnerable to subjectiveerrors [1]. Eliminating user interaction will also provide the means for real-time dosimetry, where planning the seed placement can be corrected in theoperating room during the brachytherapy procedure to account for anatom-ical changes due to patient positioning, internal organ fillings and edema aswell as seed placement error.Definition of PTV is important because it is used as a destination whereseeds are usually distributed. PTV is defined to account for intra-operativeerrors by anisotropic dilation of the CTV. Although the dilation param-eters are generally suggested in clinical guidelines, a radiation oncologistnormally modifies the contours in different anatomical regions to accountfor some implicit characteristics of the patient review chart, and his/herown expertise. The standard-of-care at the VCC is to use a manual or semi-automatic approach in contouring the target volumes, which highly dependson the expertise of a radiation oncologist.The need for CTV and PTV delineation is not limited to the preplan-ning application. Ideally, segmentation correction or plan refinement justprior to or during implantation facilitates the treatment objective: Prostateshape and size varies due to bladder and rectum fillings, cancer evolutionrate and characteristics, patient conditions, the effect of general anesthesia,prostate deformation and the development of edema during the procedure.Hence, there is a high demand for real-time target anatomy segmentation inprostate brachytherapy for clinical applications such as intraoperative pre-planning, where the plan is developed in the operating room immediately61.2. Uncertainties in LDR Brachytherapyprior to implantation, or intraoperative planning in which plan is correctedas the needles are inserted. According to the ABS recommendations forprostate dosimetry, “Ideally, one should strive for on-line, real-time intra-operative dosimetry to allow for adjustment in seed placement to achieve theintended dose” [10]. However, the acknowledged limitations of the currenttarget volume delineation process make intraoperative real-time segmenta-tion challenging. A fast and automatic delineation algorithm can potentiallylead to a more efficient and less user-dependent target segmentation, whichin turn will produce more consistent dosimetric treatment plans while facil-itating potential solutions for real-time dosimetry.Inevitably, operating room (OR) related user variability introduces asignificant amount of uncertainty to the whole treatment procedure, sincevolume study and preplanning take place weeks before the implantation.Deviation between planned and delivered dose mainly rises from three majorsources:• Prostate visibility must be the same as the volume study day. There-fore, patient, brachytherapy grid template, and the probe setup needto be replicated from the day TRUS volume is acquired. However,internal organ condition is not necessarily the same. Pelvic inter-nal muscle tension changes when patient is under general anaesthesia.Moreover, deformation and movement forced on the anatomy by theTRUS probe is not reproducible completely.• There are some technical challenges associated with delivery of theplan due to needle bending, friction, and deformation of the anatomyby the force applied through needle insertion.• The prostate boundary is not necessarily the same as the volume studyday because of the organ boundary changes. These changes are dueto cancer evolution or some medications taken meanwhile. Bladderand rectum fillings may also affect the anatomy on the operation day,although patients are instructed to empty both.Given the aforementioned sources of variability, brachytherapy treatmentoutcome and morbidity can benefit from optimization approaches that re-duce or eliminate sources of uncertainty from parts or the whole procedure.For instance, automation of the delineation processes eliminates user inter-action, however it needs to be capable of incorporating implicit knowledgeof radiation oncologists as well as clinical guidelines. A solution to such op-timization is to model the brachytherapy procedure as an integrated system71.3. Backgroundof interacting processes. This model can be used to optimally predict ele-ments of each individual procedure. The proposed model can be generatedfrom previous successful treatment records and applied to future records.1.3 BackgroundReview of the prior art in optimization (automation) of the prostate brachyther-apy procedure can be divided into two sections: 1) segmentation of the targetvolume in TRUS images; and 2) automation of the seed planning inside thePTV.1.3.1 Target Volume SegmentationAlthough there has been a lot of work on segmentation of the prostate inTRUS images, they are mostly not directly applicable to segmentation ofthe CTV and the PTV. This is further discussed in Chapter 4. To the bestof our knowledge, there is no work on automatic delineation of the PTVfrom TRUS images. Hence, in this section we review the closest works, i.e.prostate segmentation for TRUS images.Several groups have proposed semi-automatic and automatic techniquesto alleviate the manual prostate segmentation process in TRUS images. Ithas been shown that pure texture features of the TRUS images can aidclassification of the pixels and delineation of the prostate boundary [14,15]. Active contours and snakes have also been utilized by several groupsto delineate the prostate boundary in 2-D [16, 17] and 3-D [18–22] TRUSimages.Another group of algorithms take advantage of a priori knowledge ofthe prostate shape and TRUS image statistics to improve segmentation re-sults. For example, geometrical model-based approaches using super-ellipseshave been recognized as a suitable modelling approach for determining theprostate boundary [1, 23, 24]. One of the key advantages of the super-ellipses-based segmentation approaches is the strong shape regularizationthat compensates for the relatively low resolution and sparse characteris-tics of the brachytherapy TRUS volumes. In addition, this approach meetsthe requirements defined in the BCCA treatment protocol about symmet-ric CTVs. However, the proposed approaches are semi-automatic, hence,vulnerable to subjective errors and initial conditions. Although fusion ofdifferent imaging modalities has been proposed to automate the super-ellipses-based segmentation pipeline [25], it is still far from being integrated81.3. Backgroundwith current standards of care. Moreover, prior information of TRUS im-ages, such as texture features and prostate shape variability in the formof atlas or Statistical Shape Models (SSM) [26–31] have been used to im-prove robustness of segmentation. For example, Kernel Support Vector Ma-chines (KSVM) [32, 33], Active Appearance Models (AAM) [34] and levelsets [35, 36] are some of different statistical modeling approaches proposed toincorporate shape and intensity priors into a robust to noise segmentationalgorithm. In an atlas-based segmentation approach, a statistical model,i.e. atlas, is usually transformed into a common coordinate system with thetarget image. Subsequently, prior knowledge is propagated from the atlastowards the target image. Multi-atlas-based methods are based on intensityimage registration and label fusion techniques and generally have not beenused in the context of TRUS images due to the poor ultrasound image reg-istration quality. This problem also affects process of selection and fusionof atlases that keeps it challenging for the TRUS images. Therefore, singleor multiple atlas-based segmentation approaches are proposed mainly forsegmentation of the prostate in CT and MR images [37–41].Except [1, 24–26], the above mentioned TRUS segmentation methods arenot directly applicable to prostate brachytherapy without significant changesto the underlying methods or clinical workflow, because: 1) The delineatedCTV in clinical images does not follow the anatomical prostate boundaryeverywhere; 2) TRUS volumes are sparsely acquired with a slice distanceof 5 mm. This sparse nature of brachytherapy TRUS images in additionto the intrinsic artifacts of the ultrasound images, makes the segmentationprocess even more challenging. Mahdavi et al. [25] combined informationfrom another modality, i.e. elastography maps, to aid the segmentationprocess, however it requires the introduction of additional hardware or theneed for special ultrasound machines for elasticity imaging. Ghose et al. [26]proposed a solution for automatic prostate segmentation that appears to beapplicable for prostate brachytherapy. However, the approach is based onthree independent 2-D active shape and appearance models that are gener-ated only for central, apex and base zones; hence, segmentation is inherently2-D and does not guarantee a smooth CTV delineation that considers the3-D shape of the prostate.A more general approach using atlases is the multi-atlas segmentationmethod that is primarily used for segmentation of MR brain images [42–44],and more recently for CT and MR images of prostate [37, 38, 41]. Multi-atlas-based methods are based on intensity image registration and labelfusion techniques and generally have not been used in the context of TRUSimages due to the poor ultrasound image registration quality. This problem91.4. Proposed Solutionalso affects process of selection and fusion of atlases that keeps it challengingfor the TRUS images.1.3.2 Seed Arrangement PlanningThe main objective in seed planning is to target the PTV with minimumprescribed dose (mPD). Seed arrangements can be planned either prior to,right before or during the implantation procedure. The BCCA program isbased on prostate brachytherapy preplanning, which requires pre-operativeplanning of seed arrangements, usually weeks before the implantation pro-cedure.Over the last two decades, several solutions have been proposed to au-tomate seed planning [45–52]. A common approach for this automation isto define an inverse problem based on dose constraints applied to parame-ters that characterize features such as target coverage, dose uniformity, anddose to organs at risk. In this regard, various optimization techniques areemployed including simulated annealing [47], genetic algorithms [50, 52] andbranch-and-bound solution for mixed-integer model [45, 48, 51, 53]. The ini-tialization for these techniques is normally with a random seed arrangement.Considering the very large size of the search space for optimal seed config-uration, these solutions are highly subject to local minima and sensitive toinitial seed configurations.1.4 Proposed SolutionA typical brachytherapy preplanning procedure can be subdivided into dif-ferent processes as shown in Figure 1.5. Obviously, uncertainty propagatesthroughout the cascading processes and would affect quality of the treat-ment procedure, i.e. increasing rate of morbidity. Hence, current prostatebrachytherapy procedure can substantially benefit from reducing the needfor user interaction by automation of parts or the whole planning processes.In this thesis, we aim to facilitate pre-operative processes, where the currentprocedure is performed manually or semi-automatically.Current preplanning procedure is sequential and takes place in a geo-metrical extent of a patient’s TRUS volume. Hereafter, we refer to inputsand outputs of processes, i.e. TRUS intensity volume, the CTV and thePTV, as volumetric information. Since each process depends on the pre-ceding processes and observed volumetric information, we propose to fuseall available volumetric information and learn their joint relation to devise101.4. Proposed SolutionTRUS Volume Clinical Target (CTV) Planning Target (PTV) Preplan Guidelines Volume Study Figure 1.5: Current prostate brachytherapy preplanning requires user inter-vention as shown by circles.an automation solution. Our approach is to combine and capture the inter-relation between volumetric information generated from processes in theform of mathematical and statistical models.Architecture and parameters of such models are tuned and optimizedusing previous treatment knowledge. We take advantage of the previous ex-cellent outcome treatment records at the VCC in our analysis. We proposeto represent joint relation between observations of the volumetric informa-tion by joint statistical analysis of these treatment records. To this end,we investigate some statistical modelling approaches such as independentcomponent and sparsity analyses, all aiming to manifest the joint relationbetween volumetric information in a limited set of complex joint patterns.1.4.1 ObjectiveThe global objective of the thesis is to reduce the need for user interaction inplanning prostate brachytherapy treatment. To this end, we investigate au-tomation approaches in preplanning, namely clinical/planning target volumedelineation in TRUS images and seed arrangement. We investigate a fusionsystem that automatically delineates the CTV with accuracy comparable tomanual or semi-automatic methods currently used in clinics. Furthermore,we propose a fusion framework that extracts joint patterns between volu-metric information representations of the preplanning elements, i.e. CTVand PTV contours. In this thesis, we investigate the feasibility of automaticseed planning by combining the joint learning-based approach with a novel111.5. Materialsin-house optimization algorithm. Throughout this thesis, the joint analysisframework is generated and evaluated on a dataset of previously deliveredsuccessful treatment records obtained from the VCC.1.4.2 ContributionsThe present thesis is an attempt to develop a fusion-based framework to au-tomate prostate brachytherapy preplanning, according to the BCCA treat-ment program. The proposed model encapsulates variability associated withdifferent elements of determining CTV and PTV as well as seed arrange-ment using existing prior knowledge in a dataset of treatment records. Inthe course of achieving this objective, the following contributions were made:• Proposing a multi-atlas-based segmentation approach for automaticdelineation of the CTV from TRUS images.• Proposing an approach for fast and reliable estimation of the CTV,the PTV and the mPD contours using joint independent componentanalysis (jICA).• Proposing a framework for simultaneous estimation of the CTV andthe PTV using joint sparse analysis.• Proposing a seed arrangement optimizer to enforce the BCCA prostatebrachytherapy preplanning clinical guidelines.• Proposing a novel learning-based seed planning framework.1.5 MaterialsMaterials used in this dissertation are collected from the VCC under theapproval of the University of British Columbia and BCCA Research EthicsBoard, certificate number H13-01983. The target population for this studyis composed of male patients who have undergone brachytherapy treat-ment. This research requires a retrospective dataset of brachytherapy pa-tient records which includes the transrectal ultrasound images of the prostatebefore implantation of seeds, corresponding clinical target contours, plan-ning target contours, seed placement plan and plan dosimetry parameterssuch as seed activation unit, radial and geometrical dose distribution func-tions. We looked at the anonymized and de-identified data records of pa-tients who underwent treatment between May 1, 2006 and May 1, 2012. Intotal, 2,087 patient charts were reviewed.121.5. Materials(a)(b) (c)Figure 1.6: Brachytherapy grid was superimposed over the original ultra-sound images. a) A sample case captured from a video monitoring deviceexcluded from this study due to lack of an ultrasound recovery solution;and b) a sample image collected from BK machine before ultrasound imagerecovery; and c) after ultrasound image recovery.Patient brachytherapy charts had been managed and archived usingVariSeedTMLDR Treatment Planning System. All ultrasound images werealtered by the ultrasound machine software through superimposing a brachyther-apy grid and annotations. As a result, a clean up process was necessarythrough image processing software to recover the original ultrasound data.This was a challenging task, since majority of the cases in the original co-hort were captured from a video monitoring screen (some analog), hence,131.6. Thesis Outlinehad substantially degraded image quality. We focused on a smaller cohortof 590 cases where a BK ultrasound machine was used for imaging. Imagesfrom BK machines had sufficient quality to develop the recovery algorithm,and further, they have been recently widely used all across centres in BritishColumbia after 2010. For the BK machine images, a consistent procedurewas implemented to remove superimposed information and recover the orig-inal ultrasound images with minimum loss. Figure 1.6 shows some of thesample images from the cohort of this study before and after ultrasoundimage recovery.1.6 Thesis OutlineThe rest of this thesis is subdivided into four chapters as outlined below:Chapter 2: Atlas-based Fusion for Brachytherapy PreplanningIn this chapter, we introduce a multi-atlas fusion framework to automati-cally delineate the clinical target volume in ultrasound images. A dataset ofa priori segmented ultrasound images, i.e. atlases, is registered to a targetTRUS image. We introduce a pairwise atlas agreement factor that com-bines an image-similarity metric and similarity between a priori segmentedcontours. This factor is used in an atlas selection algorithm to prune thedataset before combining the atlas contours to produce a consensus seg-mentation. The proposed method produces segmentation results that arewithin the range of observer variability when compared to a semi-automaticsegmentation technique that is routinely used at the VCC.Chapter 3: ICA-based Fusion for Brachytherapy PreplanningIn search of a faster technique for fusion and estimation, in this chapter weintroduce a joint ICA-based estimator for two different applications:• Delineation of the CTV: Real-time dosimetry and intra-operative plancorrection requires a fast CTV segmentation algorithm. In this chap-ter, we propose a computationally inexpensive and fully automaticsegmentation approach that takes advantage of previously segmentedimages to form a joint space of images and their segmentations. Weutilize joint independent component analysis method to generate amodel which is further employed to produce a probability map of thetarget segmentation. We show that the proposed approach is fast with141.6. Thesis Outlinecomparable accuracy and precision to those found in previous studieson TRUS segmentation.• Delineation of the PTV and the mPD simultaneously: We introducea jICA-based model that enables joint determination of PTV and theminimum prescribed isodose (mPD) map by capturing the correla-tion between different volumetric information elements consisting oftransrectal ultrasound (TRUS) volumes, PTV and isodose contours.Taking advantage of the same jICA technique, we obtain a set of jointcomponents that optimally describe such correlation. We also performa component stability analysis to generate a model with stable param-eters that predicts the PTV and isodose contours solely based on anew patient TRUS volume.Chapter 4: Sparsity-based Fusion for Brachytherapy PreplanningIndependent component analysis can be seen as a special case of a sparsemodeling. Hence, in this chapter, we extend the ICA-based framework bytargeting simultaneous delineation of the preplanning elements, i.e. theCTV and the PTV using sparse analyses. In this chapter, we aim to re-duce the segmentation variability and planning time by proposing an effi-cient learning-based multi-label segmentation framework. We incorporatea sparse representation approach in our methodology to learn a dictionaryof sparse joint elements consisting of images, and clinical and planning tar-get volume segmentation. The generated dictionary inherently captures therelationships among elements, which also incorporates the institutional clin-ical guidelines. We show simultaneous segmentation results compared withthe currently used clinical algorithm for both target volumes.Chapter 5: Automatic Seed Plan EstimationIn this chapter, we aim to reduce the preplanning variability by automat-ing the seed arrangement process. We propose a novel framework whichuses a retrospective treatment dataset to extract common radioactive seedpatterns. The framework captures the inter-relation between the treatmentvolume delineation and seed arrangements through a joint sparse represen-tation of retrospective data. This representation is used to estimate aninitial seed arrangement for a new treatment volume, followed by a noveloptimization process which captures the clinical guidelines, to fine-tune theseed arrangement.151.6. Thesis OutlineChapter 6: Conclusion and Future WorkThis chapter includes a short summary of the thesis followed by suggestionsfor future work.16Chapter 2Atlas-based Fusion forBrachytherapy Preplanning2.1 IntroductionIn low-dose-rate brachytherapy, reliable CTV delineation in TRUS imagesis the foundation of generating appropriate seed plan to deliver sufficientradiation dose to the cancerous tissue. Current standard-of-care at the VCCrequires user intervention in determination of the CTV contours. A state-of-the-art algorithm proposed by Mahdavi et al. [1] has been integratedinto the current workflow that initiates the contours in 2-D planes. Thesecontours are further manipulated based on expert’s knowledge to ensurecancerous target is covered. This chapter is an effort to automate the CTV(approximately the prostate boundary) delineation. The main objective isto take advantage of previously contoured TRUS volumes and the state-of-the-art atlas-based segmentation approaches to facilitate prostate/CTVsegmentation process.In this chapter, we introduce an algorithm for 3-D segmentation of theCTV based on multiple atlases that include pairs of images and their seg-mented labels. The general approach of the multi-atlas segmentation (MAS)method is to transform atlases in an existing dataset to the coordinates ofa target image. Subsequently, a consensus segmentation of the target imageis produced by fusion of the atlas labels. This approach has been previouslyused for segmentation of brain MR images [42–44], and more recently forCT and MR images of the prostate [37, 38, 41, 55]. The application of MASfor TRUS segmentation is not straightforward, since it relies on robust reg-istration of images, which is a challenging task in TRUS data. To improvethe robustness to registration inaccuracy, several groups have proposed so-lutions to identify optimal weights for fusion, mainly by applying the fusionThis chapter is adapted from [54]: S. Nouranian, S. S. Mahdavi, I. Spadinger, W.Morris, S. Salcudean, and P. Abolmaesumi, AMulti-Atlas-Based Segmentation Frameworkfor Prostate Brachytherapy, IEEE Trans. Med. Imaging, vol. 34, no. 4, pp. 950-961,2015.172.2. Methodsin a neighbourhood around each target voxel, and by introducing intensitypriors to the fusion algorithms, e.g. locally weighted majority voting [56],local maximum a posteriori STAPLE [57] and local logarithmic opinion poolSTAPLE [58]. Our observation from both intensity images and their corre-sponding label maps shows a weak correlation between intensity similaritymetrics and label overlap measures in the context of sparse TRUS images.Hence, some recently introduced methods that incorporate intensity basedpriors in fusion of atlases may not perform as well in the context of ourstudy. An alternative approach proposed by Langerak et al. [59] uses apreregistration atlas selection for multi-atlas segmentation. The approachmaintains/improves the quality of input atlases based on a heuristic selec-tion of atlases prior to the registration process. The authors observe slightlybut not statistically significantly less accuracy for the target segmentation;however, the computational time was significantly reduced. In this chapter,we propose an alternative approach, where we perform atlas selection priorto and after the registration: 1) we use a target specific approach for atlasselection prior to the registration and fusion steps; 2) we introduce an at-las pruning technique to incorporate shape deformation agreement amongatlases along with the registration performance. We incorporate a pairwiseatlas agreement factor to select an appropriate subset of atlases for fusion,and subsequently produce the consensus segmentation. To satisfy CTVrequirements according to the treatment protocol at our institution, con-sensus segmentation generated by the fusion process is smoothed and madesymmetric by fitting tapered and warped ellipses. The method is evalu-ated on a clinical dataset of 280 prostate TRUS volumes. We compare ourproposed configuration for the framework against the manually segmentedgold-standard, and show it performs within a range of variability same as astate-of-the-art semi-automatic segmentation approach [23] that is part ofthe standard-of-care at the Vancouver Cancer Centre (VCC). The initial re-sults of our approach was reported in [60, 61]; here, we provide further detailsof the approach, and validate it with a much larger dataset. Moreover, wepropose to improve the robustness of the approach by introducing the atlasagreement factor. We also compare our proposed multi-atlas segmentationframework with local and intensity-driven label fusion techniques.2.2 MethodsFigure 2.1 shows a schematic block diagram of the proposed MAS frame-work. In this framework, the term atlas refers to an intensity image and182.2. Methodsits corresponding morphometric properties, such as its binary segmentation,which we hereafter refer to as the label image. Generally, MAS approachesassume a strong presence of corresponding structures in the intensity andthe label image of an atlas. They use an image registration approach toalign a group of atlases with a target image, so that the corresponding la-bels of the atlases are propagated to the coordinates of the target image. Aprobabilistic map of the consensus segmentation is achieved by fusing thetransformed labels of the atlases.The key assumption in MAS is that a perfect registration between anatlas and the target image leads to a perfect labeling of the target image bypropagating the same transformations from the intensity image to the labelimage. Where there is imperfect registration, fusion of the information frommultiple atlases reduces the sensitivity of the consensus segmentation to theregistration quality and the error associated with each individual atlas label.In the context of TRUS images obtained for brachytherapy, registrationis a particularly challenging problem, which could lead to relatively poorand spatially inconsistent results, due to the wide range of variability inshape and size of the prostate. Furthermore, fusion of atlas labels based onintensity criteria is highly affected by these specific characteristics of TRUSimages which results in an inaccurate alignment of the the atlases and hencepoor fusion performance. This implies a requirement to prune atlases beforeand after registration prior to the fusion step.Hereafter, we refer to the target volume with V , and each atlas with pair{Ii, Si}, i ∈ {1, 2, ...,M}, which represents the intensity volume Ii and thecorresponding segmentation binary volume Si : Rd → {0, 1} for the ith atlasfrom total number of M atlases in a d dimensional space.2.2.1 Target-specific Atlas SelectionIn the first step, we reduce the size of the atlas dataset using a target-specific pruning approach, where a smaller set of atlases that are similar tothe target image are selected. We follow an iterative procedure to generate amean atlas [39]. This method does not require inverting nonrigid coordinatetransformations. Subsequently, we utilize the mean atlas to bring the targetimage, along with all atlas intensity images, to a common coordinate sys-tem. Afterwards, a similarity metric between the transformed target imageand each individually transformed atlas is computed to generate a list ofcandidate atlases.By computing the image similarity metric between each pair of trans-formed atlases, {I˜ni , S˜ni }, and the transformed target image, V˜ , a fixed num-192.2. Methods... Query Image Atlas i Similarity-based Selection Mean Atlas M Registration Registration Atlas i N < M Atlas 1 Registration Deformation Atlas 2 Registration Deformation Atlas N Registration Deformation Target Segmentation Atlas Selection Label Fusion Tapered Ellipse Fit Target-Specific Atlas Selection: Filtered Multi-atlas Fusion: Clinical Target Volume Estimation: ... Atlas i M Mean Atlas Generation Figure 2.1: Block diagram of the proposed MAS-based TRUS segmentationframework.ber of N atlases expected to have a higher resemblance to the target volumeis selected. Notably, this process does not add a significant computationalburden to the segmentation pipeline as the mean atlas, I˜n, and all trans-formed atlases are computed, once, prior to the segmentation process.202.2. Methods2.2.2 Filtered Multi-atlas FusionThe second block of the proposed MAS framework consists of three cascadeprocesses; namely, registration, atlas filtration and label fusion. For theregistration, we use the diffeomorphic Demons registration algorithm [62–64], which has been shown to perform well in the context of ultrasoundimages [65, 66].For the atlas-selection, a general consensus in the MAS literature [42,44, 67, 68] is to use the image similarity metric as a measure to rank atlases.While this methodology may be effective in CT and MR images, in TRUS,due to speckles, shadowing and low signal-to-noise ratio, the sensitivity ofthe image similarity metric to anatomical variability in terms of shapes andsizes of the prostate is relatively low. To alleviate this issue, we propose anew similarity metric that jointly encodes not only the image similarity inthe intensity domain, but also the prostate shape and size similarity in thelabel domain. We define a pair atlas agreement factor, fi,j , between atlas iand j as the harmonic mean of the intensity and label distances:fi,j =2dIi,jdSi,jdIi,j + dSi,j, (2.1)where the distance between intensity of the atlases is represented by Sumof Squared Differences (SSD), denoted by dIi,j . Also, the volume error, per-centage of non-overlapping volume, is used for shape dissimilarity dSi,j . Foreach atlas, {I˜i, S˜i}, we define the pairwise atlas agreement vector Fi by cal-culating the agreement factors fi,j , j ∈ {1, ..., N} to all other transformedatlases. We follow a consistent indexing of atlases in all calculations:Fi = (fi,1, fi,2, ..., fi,N ) . (2.2)The feature vectors Fi, i ∈ {1, ..., N} are divided into K clusters G ∈{G1, ..., GK} using the k -means clustering algorithm. Intuitively, membersof each cluster share a similar distance pattern with the other members.We utilize the cosine similarity between pairwise atlas agreement vectors asthe distance metric for the clustering algorithm. The k -means algorithmminimizes the within cluster distance between members Fi and the centroidof each cluster C1, ...,CK :arg minGK∑i=1∑Fj∈Gi(1− Fj .Ck‖Fj‖‖Ck‖). (2.3)212.2. MethodsTo minimize the effect of initialization on the clustering algorithm, cen-troids of the clusters are randomly initiated many times and the best clus-tering result with the minimum summation of the distance to the centroidsis selected. Subsequently, clusters are ranked based on their average imageintensity similarity to the target image. We define the ranking parameterrGi for cluster Gi as the mean SSD between members of cluster Gi and thetarget image:rGi =1|Gi|∑I˜∈Gi∑x∈I˜(I˜(x)− V (x))2, (2.4)where x represents a voxel in the transformed atlas image I˜.Next, members of the top ranked cluster, GT , are selected for the labelfusion step. Transformations obtained from the registration step are prop-agated to the corresponding atlas labels, producing S˜i, i ∈ {1, ..., P}. Thisgroup of transformed labels are combined together to provide a consensussegmentation L̂ for the target volume V .There are several label fusion strategies introduced to determine the op-timal weight for each atlas label. A group of algorithms adapt a determinis-tic voting strategy, such as majority voting and globally weighted majorityvoting, and assign equal or adjusted weights based on a global atlas qual-ity measure (e.g., the registration performance metrics). Another group ofalgorithms include a stochastic model within a statistical fusion strategy. Si-multaneous Truth and Performance Level Estimation (STAPLE) [69] is oneof the widely used statistical fusion approaches, which has been shown toperform well for MAS-based segmentation [44, 70]. STAPLE was originallydeveloped to account for intra-rater variability by estimating an underlyingground-truth from multiple segmentations of the same object. To accountfor spatially inconsistent registration quality, locally weighted variations ofthe voting [56, 67] and local STAPLE algorithms [57, 58, 71] have been pro-posed. In our experiments, we evaluated different label fusion techniquesusing a fixed number of input atlases to produce a consensus segmentation.2.2.3 Clinical Target Volume EstimationAccording to the VCC treatment protocol [6], the CTV contours are desiredto be symmetric and smooth. The treatment protocol instructs to drawsymmetric contours with respect to the mid-sagittal plane. This is to min-imize the effect of intra-operative image reproducibility error, that causesdeviation from the pre-planned dose distribution. To accomplish this, we222.3. Materialsfollow the approach proposed by Mahdavi et al. [1] and fit a tapered el-lipse to unwarped contours. The algorithm comprises the following steps:1) obtaining the unwarping parameters based on anatomical landmarks, 2)unwarping the contours w.r.t. the unwarping parameter, 3) fitting a taperedellipse to unwarped contours, and 4) warping the fitted tapered-ellipse backto the target coordinates. The required steps are shown in Figure 2.2. Threeinitialization landmarks are identified automatically for each slice from thesegmentation result (i.e. STAPLE algorithm output) defined as 1) P1: probecenter, 2) P4: mid-posterior and 3) P2: lower posterior-lateral. The tapered-ellipse fit process is expressed as an optimization problem using the l-bfgsoptimization algorithm [72]. The explained pipeline will generate a smoothand symmetric contour with respect to the sagittal plane for all 2-D TRUSplanes.2.3 MaterialsA dataset of 280 TRUS volumes of brachytherapy patients is used for eval-uation of the proposed segmentation algorithm. All images were acquiredat the VCC, a few weeks before the treatment day. Data acquisition andplanning procedures were carried out in accordance to the local treatmentprotocol [6]. Instructions in the protocol standardize all TRUS volumesto be trimmed from prostate base to apex, while the gland is visually lo-cated in the middle of each axial B-mode image. As mentioned earlier, eachTRUS volume consists of 7 to 14 parallel equally spaced (5 mm apart) ax-ial B-mode images of the prostate which are captured using a side-firingtransrectal probe. For each B-mode image, the prostate gland is delineatedusing a software. This software is part of the standard of care, and worksbased on a state-of-the-art semi-automatic prostate segmentation methodpresented by Mahdavi et al. [1]. These contours are manually correctedafterwards by an expert clinician, and subsequently are used as the CTVreference contours for evaluation of the proposed algorithm.Images are all preprocessed to remove overlaid device labels and marks,and are down-sampled by a factor of three (136×165 pixels) to accelerate thecomputational processes. The intensity range is normalized across all atlasesw.r.t the target image by histogram matching prior to the atlas registrationprocess. Subsequently, all atlases are geometrically aligned, and if necessary,resampled (with the nearest-neighbour technique) w.r.t. the target imageextent. The extent of each individual volume is limited to the base and apexaxially, and is standardized in two other axes (sagittal and coronal) among232.3. Materialsp1p4p2(a) (b)Fitted tapered ellipsexy Original ContourTapered Ellipse Fit(c)Sample points for t−ellipse fitting Original ContourTapered Ellipse Fit(d)Figure 2.2: Process of making contours smooth and symmetric. a) initialpoint required for unwarping parameter estimation, b) unwarped image andcontour, c) tapered ellipse fit, and d) final warped tapered ellipse and theoriginal contour. All contours are generated for a sample mid-gland 2-DTRUS image.242.3. Materialsall atlases.Hyper-parameter K, which refers to the optimal number of clusters, istuned by performing validation on a subset of the dataset. Statistics of thesegmentation outcome on the validation dataset are obtained for differentvalues of K, and the optimum number of clusters is determined. Poten-tially, this parameter exposes a characteristic of the training dataset, e.g.variation in anatomy coverage, which is intuitively expected to be the samebetween validation and test datasets. We adapt a systematic sampling ofthe dataset to divide it into 60% train, 20% validation, and 20% test cases.The partitioning process is repeated 10 times in order to evaluate robustnessof the proposed segmentation algorithm. While none of the test cases areincluded in training and validation sets within a partitioning process, eachtest atlas might appear at most twice in test datasets among partitioningattempts. In each trial, the target volume is registered to the mean atlas asexplained in Section 2.2.1 in order to prune the training dataset into a finalset of 100 atlases, which are further transformed into the coordinates of thetarget image.2.3.1 Segmentation EvaluationFor every case, segmentation performance is evaluated against the clinicalgold-standard as well as the semi-automatic segmentation algorithm of Mah-davi et al. [1]. We perform statistical non-inferiority analysis using pairedt-test. We define the null hypothesis as both semi-automatic and proposedalgorithms produce results with the same mean value for the evaluationmetrics (errors). Contours of the semi-automatic algorithm are directly ob-tained from the software used as part of the standard of the care at theVCC along with the gold-standard contours. We carry out computation ofvolumetric and surface-based distance measures to evaluate the segmenta-tion performance. We use the Dice Similarity Coefficient (DSC) to quantifyshape similarity between two binary volumes, i.e. the segmentation outcomeand the gold-standard. We also denote the volume error by Verr and defineit as 100 − DSC. If we denote L̂ and L as the estimated and true binaryCTV volumes of the target image, we compute the DSC asDSC =2∣∣∣L̂ ∩ L∣∣∣∣∣∣L̂∣∣∣+ |L| × 100. (2.5)In brachytherapy, volume is a key parameter in treatment planning anddosimetry analysis, as there is a high correlation between the number of252.3. Materials Base Apex Posterior Anterior Mid-gland Lateral Figure 2.3: Prostate sectioning for segmentation evaluationseeds and needles required and the volume size. We quantize this differenceby calculating volume difference coefficient asVdiff =(∣∣∣L̂∣∣∣− |L|)|L| × 100, (2.6)which is a signed coefficient representing size difference between estimatedand the gold-standard CTV binary volumes.For the mid-gland slice contours, we represent the distance by calculatingMean Absolute Radial Distance (MAD), and Maximum Absolute RadialDistance (MAXD) in between. The mid-gland slice is expected to providethe clearest view of the prostate anatomy for the delineation process. Wealso quantify the overall surface distance by calculating the mean surface(MSD) and Hausdorff distance (HD) between two contour sets.Multi-atlas segmentation is highly dependent to the registration perfor-mance, i.e. global SSD between the registered atlases and the target images.Rohlfing [73] criticized tissue overlapping metrics and their correlation withregistration accuracy. He showed that a good global overlapping metricdoes not necessarily correspond to a good registration result. To accountfor the weakness of global metrics, we break down computation of the sim-ilarity metrics (both volumetric and surface-based) into several anatomicalregions similar to Mahdavi et al. [1]. From a clinical perspective, each of theanatomical regions are subject to a different expected accuracy and impor-tance with respect to the treatment planning protocol. These regions aredefined by partitioning the whole gland into six sectors. Axial partitioning262.4. Experiments and Resultsis made by dividing w.r.t. 0.3, 0.4 and 0.3 of the total prostate length frombase to apex. Each partition is further divided into lateral (left and right),posterior and anterior sections by two perpendicular sagittal planes passingthrough the axial line of the prostate, producing nine general anatomicalsectors as shown in Figure 2.3.2.4 Experiments and ResultsTRUS images provide a field of view that encompasses some of the adjacenttissues, including portions of the bladder and rectum. Obviously, the globalregistration metrics do not represent the registration performance related tothe region of interest. Therefore, we define a mask volume M on an areawhere there is a disagreement between transformed atlas labels:M(x) ={1, ∃(i, j), S˜i(x) 6= S˜j(x)0, otherwise(2.7)We individually obtain and use a volumetric mask M for each new tar-get image to include pixels with higher uncertainty in calculation of theregistration quality metric.All atlas and target registrations are performed in 3-D on the wholeTRUS volume. We use the Demons registration algorithm [62–64] with aGaussian kernel size of 5 mm in the axial direction.A linear correlation analysis of the intensity and label structure match-ing is illustrated in Figure 2.4. We use SSD for quantifying the intensitymatching, and DSC for binary label shape similarity. We obtain a low cor-relation coefficient of 0.29 with an arbitrary distribution of the observedpaired metrics (SSD, DSC) for all 10 experiments on partitioned dataset,which points to the necessity of an atlas pruning technique prior to and afterregistration.We choose N = 100 atlases to form the MAS framework dataset priorto the registration process. After the registration step, we apply the clus-tering on pairwise atlas agreement factor vectors and choose the highestranked cluster w.r.t the average global SSD of its members. All members ofthe winning cluster are fused using the STAPLE algorithm to produce theprobability segmentation map of the target image. Calculation of intensity-based similarity metrics is confined within the area of the generated maskfrom all transformed atlases, as explained above.272.4. Experiments and ResultsFigure 2.4: Correlation analysis between SSD (registration performance)and DSC (segmentation overlap) performed on test cases. This analysishas been achieved by running the experiments on partitioned dataset for 10times.2.4.1 Mean Atlas GenerationFigure 2.5.c shows the mid-gland slice of the average atlas generated fromthe intensity volumes of atlases in a training dataset. Blurriness of the meanatlas relates to the performance of the registration algorithm in the train-ing dataset, whereas in the current context, it is highly affected by intrinsiccharacteristics of the TRUS images. We use the Demons registration al-gorithm, and apply the mean atlas generation algorithm for 10 iterations.We found this number sufficient to converge to a consistent average volume.Figures 2.5.a and 2.5.b illustrate the mid-gland slice result of a registrationprocess between a sample TRUS volume and the mean atlas.2.4.2 Global and Sector-based Performance MeasuresTable 2.1 shows the error statistics calculated for all 10 trials, in which atotal number of 560 bootstrapped cases is included in the error calculation.We report mean and standard deviation of the error metrics to compare282.4. Experiments and ResultsFigure 2.5: For the mid-gland slice: a) a sample atlas, b) registered to themean atlas, and c) the generated mean atlas at the same slice.Table 2.1: Volumetric, contour and surface error statistics (mean±standarddeviation) calculated from 10 random dataset partitionings are shown forthe proposed automatic segmentation approach and compared with the ap-proach of Mahdavi et al. [1].Volumetric Error(%)Mid-gland ContourError (mm)Surface Error (mm)Verr Vdiff MAD MAXD MSD HDAutomatic 10.25±3.74 -2.37±13.73 1.77±0.98 4.52±2.04 1.33±0.60 5.48±1.71Semi-auto. 8.81±5.60 3.56±14.66 1.20±1.21 3.36±2.56 1.11±0.81 5.57±2.28accuracy and precision of the proposed automatic segmentation method.Since the gold-standard segmentation is based on manual initialization ofthe semi-automatic segmentation of Mahdavi et al. [1] in the mid-gland,our mid-gland contour metrics show the most significant error amongst themetrics.In Tables 2.2 and 2.3, we show evaluation metrics (i.e. DSC and MSD)obtained from 10 different random partitionings of the dataset (M = 280)into 60-20-20% training-validation-test sets per sector. To provide detailedanalysis of accuracy distribution, a volumetric measure of DSC is shownin Figure 2.6. In the posterior region, the automatic algorithm providesa higher DSC value in comparison to the semi-automatic algorithm. Incontrast, at the anterior region there is less accuracy obtained comparedto the semi-automatic algorithm. In the rest of the regions both methodsperform with similar accuracy. Noticeably, there are cases for which thesemi-automatic algorithm presents the highest DSC value, meaning therehas been no need to modify the generated contours.Two typical segmentation results are shown in Figure 2.9, each overlaidby the gold-standard and semi-automatic based contours.292.4. Experiments and ResultsPosterior Lateral Anterior Posterior Lateral Anterior Posterior Lateral Anterior50556065707580859095100DSC Base Mid ApexAutomaticSemi-autoFigure 2.6: DSC values computed for nine anatomical sectors of the prostateand for two segmentation methods: 1) the proposed automatic approach,and 2) the semi-automatic method. Statistics are computed from 10 randompartitionings of the dataset into 60-20-20% train-validation-test sets, i.e. 560test cases.New TRUS Pre-processing ~3 sec Register to Mean Atlas ~4 sec Filter Atlas Dataset ~1 sec Atlas Registration ~120 sec Atlas Pruning ~20 sec Atlas Fusion ~5 sec Post-processing ~10 sec CTV Figure 2.7: Average operating computation time measured for each individ-ual process of the proposed segmentation pipeline and its configuration withtotal number of 100 atlases.302.4. Experiments and ResultsSTAPLE Spatial STAPLE LWV SIMPLE l-MAP STAPLE l-LOP STAPLE Proposed Method5101520253035Verr (%)STAPLE Spatial STAPLE LWV SIMPLE l-MAP STAPLE l-LOP STAPLE Proposed Method-50050Vdiff (%)Figure 2.8: Comparison between different fusion algorithms on all pre-selected atlases and the proposed algorithm. Verr and Vdiff represent shapeand size dissimilarity between the estimation result and the gold-standardfor all test cases in a random dataset partitioning.312.4. Experiments and ResultsTable 2.2: Grand mean and pooled standard deviation of Dice (DSC) cal-culated for 10 different trials of randomly dividing data into 60-20-20%training-validation-test sets.AutomaticDSC % Base Mid Apex TotalAnt. 88.40±8.50 90.24±7.84 85.78±10.7889.75±3.75Lat. 91.26±4.62 91.42±4.56 87.54±6.04Post. 91.48±4.87 92.12±4.75 87.33±8.13Semi-automaticDSC % Base Mid Apex TotalAnt. 92.22±8.69 93.01±7.20 87.70±12.7491.19±5.60Lat. 92.96±6.02 93.35±5.57 87.47±10.17Post. 90.42±7.63 94.35±4.75 88.81±11.942.4.3 Computational CostStandard-of-care at the VCC, where the TRUS data for this study werecollected, is to conduct the planning process a few weeks before the sched-uled implantation date, a time interval during which the treatment plan isgenerated and the seeds and needles ordered from the manufacturer. Cur-rently, 3 ∼ 4 min of the planning procedure is dedicated to the segmentationprocess [1], and no time constraints are applied to the segmentation processas it is a very early step in the planning procedure. However, because itis amongst the first of several steps in the planning pipeline, it should beperformed as quickly and efficiently as possible.We measure execution time of the proposed automatic segmentation ap-proach to ensure there is no burden added to the standard of the care atthe cost of automating the segmentation procedure. Computation time ismeasured on a standard PC (Intel Core i7, 2.93GHz, 8GB RAM). All themain functions are implemented in C++ and MATLABR©.The mean atlas generation process takes place only once. The mean atlasand all atlases transformed to the coordinates of the mean atlas are storedfor the subsequent steps. Fig 2.7 provides an overview of the computationtime of individual processes of the proposed algorithm. Reported numbersare for the non-optimized code employed in two different programming in-terfaces (i.e. C++ and MATLABR©). The proposed method runs under 3minutes according to the maximum execution time measured for the evalu-ation dataset. Therefore, the automatic approach is fast enough that allowsthe physicians to perform more tasks such as contour modification while it322.4. Experiments and ResultsTable 2.3: Grand mean and pooled standard deviation of mean surfacedistance (MSD) calculated for 10 different trials of randomly dividing datainto 60-20-20% training-validation-test sets.AutomaticMSD mm Base Mid Apex TotalAnt. 0.96±0.78 0.82±0.70 0.65±0.661.33±0.60Lat. 0.92±0.56 0.96±0.56 0.69±0.47Post. 0.63±0.38 0.56±0.34 0.53±0.44Semi-automaticMSD mm Base Mid Apex TotalAnt. 0.58±0.69 0.63±0.60 0.51±0.641.11±0.81Lat. 0.68±0.66 0.82±0.71 0.58±0.55Post. 0.59±0.52 0.37±0.30 0.45±0.60generates contours for a patient.2.4.4 Comparison of Different Fusion MethodsFigure 2.8 shows volume and shape distance metrics obtained by the com-pared fusion algorithms. A total of 100 atlases, preselected for the multi-atlas registration and fusion, are processed by: 1) the conventional STAPLEalgorithm [69]; 2) spatial STAPLE [71]; 3) Locally Weighted Majority Vot-ing [56, 67] (LWMV); 4) Selective and Iterative Method for PerformanceLevel Estimation (SIMPLE) [68]; 5) local Maximum A Posteriori STA-PLE (l-MAP STAPLE) [57]; and 6) local Logarithmic Opinion Pool STA-PLE [58] (l-LOP STAPLE). We have used implementation of the local LOPand the local MAP STAPLE from http://crl.med.harvard.edu/, and theSIMPLE and spatial STAPLE from http://www.nitrc.org/projects/masi-fusion/. The parameters used were those recommended by the authors ofthe original algorithms. In local variation of the algorithms, we limit neigh-bourhood window size to 1.4 mm in 2-D axial planes. Statistical significanceanalysis of the measured Verr using paired t-test shows significant improve-ment by the proposed method when compared against the evaluated fusionalgorithms (p-value< 0.01). Among the other fusion approaches, local LOPSTAPLE provides the most accurate segmentation results both in terms ofshape and size, as measured by Verr and Vdiff .332.4. Experiments and ResultsBase Mid-gland Apex (a)Base Mid-gland Apex (b)Figure 2.9: Two examples of the automatic segmentation result comparedto the semi-automatic approach of Mahdavi et al. [1]. The gold standardcontour is shown by dashed line.342.5. Discussion and Conclusion2.5 Discussion and ConclusionWe presented an automatic 3-D segmentation algorithm for delineation ofthe CTV in prostate brachytherapy. The algorithm uses a non-parametricdeformable registration method to align a set of a priori atlas volumes with atarget volume. The CTVs delineated for all atlas volumes are then combinedin a STAPLE framework to provide a probability map of the target imageCTV. We proposed post-processing steps by adapting a tapered-ellipse fitalgorithm to obtain the consensus segmentation of the CTV. The proposedalgorithm can potentially reduce the inter- and intra-observer variability ofthe CTV contouring by eliminating the need for user intervention. How-ever, modification of the automatically generated contours by clinicians isinevitable for seamless translation of the approach to clinic.To reduce the computational cost and maintain the efficiency of thealgorithm, we take advantage of a forward mean atlas generation processprior to the registration of atlases. We filter a large dataset of atlases into asmaller set of similar atlases using the generated mean atlas, hence reducingthe chance of participation for those atlases with lower similarity betweentheir intensity volumes and the target volume.Our observation from the registration performance concludes that par-ticipation of all atlases in the fusion procedure decreases performance andaccuracy of the segmentation; namely, reduces the mean DSC and raises thestandard deviation. One reason can be the fact that registration accuracyis not equally spread over different anatomical regions within the region ofinterest. Moreover, there is no strong evidence behind the assumption of thelower SSD in intensity volumes, the higher DSC between binary volumes isexpected. For the ultrasound modality, deformation fields obtained fromregistration algorithms are highly affected by the presence of speckles andother artifacts that appear as registration noise in the boundary region ofthe prostate. Therefore, propagation of the deformations from the intensitydomain to the label’s binary domain is subject to some distortions. Even-tually, filtering the registered atlases instead of including all atlases reducesthe registration error.In a major departure from earlier prostate segmentation approaches, theproposed pipeline is applicable to the sparse ultrasound volumes collected inbrachytherapy. Intuitively, we expect a better registration quality for densevolumes, because more anatomical features contribute to the registrationand fusion processes. However, this will increase computational cost of thesegmentation pipeline.A non-optimized implementation of our proposed algorithm runs under352.5. Discussion and Conclusion4 minutes in a standard PC, which is comparable with the current segmen-tation speed in a typical brachytherapy volume study process [1]. However,individual atlas registration to the target volume is an independent process,and a parallel implementation on multiple CPU cores and possibly GPU canreduce the computation time further to the order of seconds. Therefore, ourmethod has the potential to handle the requirement of real-time dosimetryfor brachytherapy with intra-operative planning.2.5.1 Comparison with Other MethodsA comparison analysis with other advanced fusion algorithms shown in Fig-ure 2.8, indicates that the low correlation between structures in intensityand label domains can mislead the local or intensity-driven weight assign-ment performance, and hence may undermine the fusion results. This ob-servation implies that selecting the appropriate atlases prior to fusion cansignificantly improve the accuracy of consensus segmentation in the contextof TRUS images. Our framework can potentially benefit from other atlasselection methods, such as the approach recently proposed by Asman etal. [74] which has been evaluated on MR images of cervical vertebrae.We retain clinically comparable results in terms of accuracy and preci-sion against the gold-standard CTVs on a dataset of sparse TRUS volumesobtained from 280 patients. The inter- and intra-observer variability ofthe prostate segmentation in the context of brachytherapy has been doc-umented before [1]. Specifically, the inter-observer variability of manuallysegmented contours has been reported to be on the order of 7.25 ± 0.39%and 6.64 ± 2.36% for Verr (i.e 100 − DSC) and volume difference (Ddiff ),respectively. Moreover, the intra-observer variability of manually segmentedcontours has been reported to be on the order of 5.95 ± 1.59%. Our pro-posed fully-automatic method improves robustness and consistency of thesegmentation, and produces results within the accepted range of variability.We also compare segmentation results to some reported results in theliterature; however, direct comparison is not possible. Tutar et al. [75] re-ported volume overlap of 82.8±6.2% (corresponding to DSC value of 90.59%)as disagreement between three different experts in manual boundary delin-eation of 30 post-implant TRUS prostate volumes. Akbari [32] and Yang [33]obtained DSC values of 88.1 ± 1.44% and 90.7 ± 2.5%, respectively, for asmall dataset of five dense TRUS volumes. Heimann et al. [34] reportedoverlap error of between 16.7± 5.2% to 17.6± 6.8% (corresponding to DSCvalue of 90.35% and 90.89%) for successful segmentation of 25-33 out of 35dense TRUS volumes. Gong et al. [23] reported inter-observer variability362.5. Discussion and Conclusion50 60 70 80 90 100020406080100DSC(%)test cases (%) 85AutomaticSemi-automaticFigure 2.10: Performance comparison between the proposed automatic ap-proach and the semi-automatic method presented as percentage of test casessegmented with respect to DSC values. The proposed method achievesgreater success rate for DSC values less than 86% compared to the semi-automatic method.of 1.82± 1.44 mm for the contour distance in all planes obtained from fiveexperts and out of 125 images of 16 patients.Figure 2.10 provides a comparison chart for segmentation success rate, interms of DSC with the gold-standard, on a 10 random dataset partitioning.Herein, we define the accuracy threshold of DSC as a parameter to computethe success rate, i.e. percentage of results with greater DSC value. It canbe observed that for the accuracy threshold of 85% the proposed MAS ap-proach is successful for 90% of the test volumes (i.e. 560 cases), which iscomparable to the success rate of the semi-automatic algorithm [1] whichis 87%. However, the semi-automatic approach outperforms in terms of thesuccess rate for higher segmentation accuracy threshold values.Detailed evaluations on nine anatomical sections shows maximum sur-face error in the posterior-apex section (see Figure2.6). Although averageDSC obtained from the semi-automatic algorithm is higher than the pro-posed MAS framework, less standard deviation is seen in most of the re-gions. Highest DSC values are obtained in the mid-gland section in bothalgorithms.2.5.2 Dataset DependencyDefinition of CTV and subsequently the Planning Target Volume (PTV)is a subjective process which varies in between brachytherapy clinics. The372.5. Discussion and Conclusionmulti-atlas segmentation approach takes advantage of existing segmenta-tions by projecting them onto the target TRUS volume. This makes theproposed MAS approach a unique method that combines intensity-basedsegmentation methods with a learning-based technique that is capable oflearning delineation style in the atlas dataset. In contrast to those segmen-tation methods which rely on boundary region detection (e.g., edge-basedmethods or active shape models), the multi-atlas approach projects the ex-isting pattern in the atlas dataset onto the target image. This feature makesthe approach a potential solution for problems similar to the brachytherapyCTV delineation, in which the clinical label does not necessarily follow thereal target anatomy boundary in all regions.A drawback of the proposed algorithm is the dependency on the dataset,which affects the accuracy of the final segmentation. However, it might bedesirable to form a generic atlas dataset that includes atlases representingdifferent prostate shapes and sizes.In summary, the proposed multi-atlas segmentation approach is a clini-cally adaptable algorithm that can be integrated with the current brachyther-apy volume study procedure. Although the segmentation execution timeis not a design criteria for the brachytherapy preplanning, the proposedmethod can be accelerated and optimized for real-time clinical diagnosticor treatment applications such as targeted biopsy and real-time dosimetry.A potential solution is to replace the Demons registration method with theGPU-accelerated version of NiftyReg [76].38Chapter 3ICA-based Fusion forBrachytherapy PreplanningIn search for a machine learning fusion algorithm that can efficiently extractjoint variations between different volumetric information representations ofthe brachytherapy preplanning elements, we investigate joint independentcomponent analysis (jICA) in this chapter. This chapter is subdivided intotwo parts aiming to utilize jICA methodology in estimation of contours; andprediction of both planning target volumes and the mPD isodose contours.3.1 Clinical Target Volume Estimation fromTRUS Images3.1.1 IntroductionIn prostate brachytherapy, usually TRUS volumes are axially sparse. Con-sidering some intrinsic characteristics of the ultrasound images and poorvisibility of the prostate near the base and the apex, CTV delineation re-mains a challenging problem. Moreover, there is a high demand for anefficient and fast prostate segmentation approach that can be used intra-operatively. This potentially allows the physicians to correct the plan w.r.t.the operating room conditions to ensure CTV is properly targeted. How-ever, most of the existing segmentation approaches are limited in terms ofthe speed and performance to fulfil intraoperative requirements includingour approach on using multiple atlases explained in Chapter 2 [54, 60, 61].Therefore, a fast segmentation algorithm with minimum interaction withthe operator facilitates potential solutions for the real-time dosimetry.For the first time in the context of prostate segmentation, we jointlyanalyze a priori knowledge of the CTV and TRUS images by generatingThis section is adapted from [77]: S. Nouranian, M. Ramezani, S. S. Mahdavi, I.Spadinger, W. J. Morris, S. E. Salcudean, and P. Abolmaesumi, Fast Prostate Segmen-tation for Brachytherapy based on Joint Fusion of Images and Labels, in Proceedings ofSPIE Medical Imaging, 2014, vol. 9036, p. 90361A1-7.393.1. Clinical Target Volume Estimation from TRUS Imagesa model using joint Independent Component Analysis (jICA). We use adataset of TRUS images and their delineated CTVs determined by expertclinicians that comprises a large variety of prostate shapes. We form a jointdata matrix including intensity images and their corresponding labels inorder to seek statistically independent basis functions, e.g. joint image-labelsources, within the existing dataset.The model generated from jICA can be potentially used in a real-timesegmentation application. The target image is projected onto the joint basisvectors in order to obtain the optimum coefficients for reconstruction. Weassume the target label can be reconstructed by a linear combination of theobtained basis vectors. Accordingly, same coefficients are employed to re-construct a probability map of the target label. The target label is estimatedfrom this probability map using a globally optimum threshold followed by acontour smoothing process. We evaluate this method against the manuallysegmented prostate contours that were used by clinicians to plan brachyther-apy procedures for 60 patients. We show clinically acceptable results thatsatisfies speed requirement for a real-time clinical application.3.1.2 MethodsModel Generation based on joint ICAIndependent Component Analysis (ICA) is a statistical approach to decom-pose a set of observations, x, into components which are maximally inde-pendent:x = aTC (3.1)where a is the mixing vector and C = [c1c2...cK ]T is the component matrix.ICA works with higher order statistics of the observation data and assumesthat hidden sources are statistically independent, non-Gaussian with lin-ear mixing process. For the context of this study, ICA is performed ona reduced dimensionality observation matrix using Principal ComponentAnalysis (PCA).Joint Independent Component Analysis (jICA) [78], is an extension ofICA used as a data fusion technique to combine information from multi-ple modalities. jICA has been successfully applied in cognitive and brainstudies [79–81]. We employ jICA to learn the correlation between intensityimages and their binary labels. We obtain a set of joint components consistof the intensity and the corresponding label maps.For each patient data, we form a joint observation vector xi that is de-scribed as 1-D representation of TRUS volume and the corresponding seg-403.1. Clinical Target Volume Estimation from TRUS Imagesmentation. A joint observation matrix X is defined as a stack of observationvectors:X = [x1x2...xN ]T where xi =[xvixli], i = 1, 2, ..., N. (3.2)Elements of the vector xvi are scalar values of a particular voxel inside theTRUS volume, and xli is a binary vector representing the label volumes.The proposed segmentation pipeline consists of two main parts, mod-elling and evaluation. The modelling process involves decomposition ofthe observation matrix X into K number of joint independent componentsC = [c1c2...cK ]T . The evaluation process uses the joint basis vectors ona test TRUS volume. The test volume is partially projected onto the sub-components cvj to obtain the mixing vector a. Same mixing vector is furtherused to combine the second part of the joint components, clj to generate asegmentation probability map. A heuristic threshold value is used to con-vert the label probability map into a binary map and its corresponding edgecontours. The generated contours are smoothed to produce the final esti-mation of the target segmentation. Figure 3.1 shows block diagram of theproposed segmentation pipeline.Label Map Reconstruction and SmoothingA label map is produced from a linear combination of the label basis vectorscj , using same coefficients obtained from the intensity image projection a.In a post-processing step, we convert the produced label map into the finalestimated binary volume using an empirical threshold value (τ = 0.4 fora normalized label map). Edge contours of the estimated binary volumeare smoothed slice by slice. Contours are represented in a polar coordinatesystem assuming their geometrical center as the origin and further trans-formed into the Euler coordinate system to perform B-spline curve fitting.Figure 3.2 shows smoothing process for a single contour of a sample result.Coordinates of the contour points are mapped to B-spline curve-fits in twovertical (x) and horizontal (y) directions separately.ExperimentsThe proposed method is evaluated on a dataset of 60 TRUS volume imagesof patients undergoing brachytherapy using leave-one-out cross validationapproach. Each TRUS volume consists of 7 to 14 transverse plane 2-D ul-trasound images axially acquired from base to apex while the prostate gland413.1. Clinical Target Volume Estimation from TRUS ImagesjICA Joint Spatially Ordered Image-Label Matrix Intensity Label Joint Decomposition Intensity Sources Label Sources Target Image Construction Post-processing Label Map Estimated Target Segmentation Model Generation (offline ) Evaluation (real-time) Projection Figure 3.1: Block diagram of the proposed segmentation approach. Modelgeneration is performed off-line and based on jICA. Real-time applicationutilizes the generated model to estimate for the target CTV delineation.is visually aligned in the middle of all ultrasound images. For each TRUSvolume, the delineated prostate gland in all 2-D ultrasound images is pro-vided by an expert radiation oncologist as the gold-standard in this study.The accuracy and precision of the algorithm is evaluated using 1) volumet-ric measures such as dice similarity coefficient (DSC) and percent volumedifference, Vdiff [1]; 2) contour distance measures such as the Mean andMaximum Absolute radial Distances (MAD and MAXD) and 3) surface dis-(a) (b) (c) (d)Figure 3.2: Contour smoothing process, (a) a sample contour, (b) x and ycoordinate coefficients, (c) B-spline curve fitting, (d) smoothed contour.423.1. Clinical Target Volume Estimation from TRUS Imagestance measures such as Hausdorff Distance (HD) and Mean Surface Distance(MSD) against the gold standard determined by the expert clinician.In order to generate the spatially ordered joint image-label matrix, allprostate image and label volumes are re-sampled towards the size of theaxially longest TRUS volume. To speed up the process, 2-D ultrasoundimages of each volume are down-sampled by a factor of three. All intensityimages are normalized using grand mean and standard deviation computedfrom dataset members. Joint image-label vectors are formed using a uniquespatial order and stacked to form the joint data matrix. Data dimensionalityis reduced to the order of 35 using PCA in order to retain over 95% ofthe information from eigenvectors. Followed by a whitening process, i.e.subtraction of the mean data, jICA decomposition is performed to obtain35 joint basis vectors.3.1.3 Results and DiscussionStatistics of the volumetric, contour and surface distance evaluation metricsare compared with the gold standard segmentation and shown in Table 3.1.A direct comparison with the reported results in the literature is not possible.However, our proposed method produces results with comparable volumetricmeasures against other TRUS image segmentation reports. The presentedresults have to be considered within the context that the gold-standard itselfis subject to inter- and intra-observer variability on top of the variability as-sociated with the segmentation method. Our gold-standard segmentationsare obtained from actual delivered treatment plans which are produced byone expert clinician. Since our comparison is with respect to manual seg-mentations of one clinician, we refer to the literature and review the clinicalsignificance of the obtained results. For segmenting the prostate bound-ary for planning brachytherapy [1], the inter- and intra-observer variabilityof manually segmented contours in terms of the volume error (1 − DSC)has been reported to be in the order of 4.65 ± 0.77% and 5.95 ± 1.59%,respectively. In another study [75], volume overlap of 82.8 ± 6.2% is mea-sured when estimating for the inter-observer variability of manual prostateboundary delineation which corresponds to 9.41% volume error. The valueof 1.82±1.44 mm [23] is also reported for the contour distance inter-observervariability in a 2-D segmentation study.Figure 3.3 shows generated contours in six slices from base to apex fortwo average segmentation outcomes based on the calculated DSC value.Volumetric (DSC) and surface distance (MSD) evaluation metric histogramsare shown in Figure 3.4 as well as TRUS images of an outlier with lowest433.1. Clinical Target Volume Estimation from TRUS ImagesTable 3.1: Comparison of various evaluation metrics between the proposedsegmentation approach and the manually segmented contours determinedby an expert clinician.Volumetric Measures (%) Contour Errors (mm) Surface Errors (mm)DSC Vdiff MAD MAXD MSD HD90.28± 3.07 −1.69± 11.85 1.84± 0.93 4.52± 1.77 1.30± 0.48 5.39± 1.46DSC and highest MSD value. Due to strong appearance of the ultrasoundartifacts, delineation of the prostate boundary is very difficult for this casein most of the slices. In addition, size of the prostate is largest among thedataset and most likely is not learned by the modelling approach whichjustifies failure of the proposed automatic segmentation approach for thiscase.The proposed segmentation procedure requires minimum computationalcost to produce the label map, since projection of the target intensity imageis formulated as a least squares problem. We measured execution time froma non-optimized code on a standard PC (Intel Core i7, 2.93 GHz, 8 GBRAM) to be 645 ± 24 ms. Our smoothing approach consists of coordinatesystem transformation and B-spline fit which adds 1158 ± 432 ms to thesegmentation pipeline. We conclude to execution time of less than twoseconds on a standard PC which is sufficient for enabling real-time dosimetryplanning.(a)(b)Figure 3.3: Contours for six equally-spaced slices (base to apex) of twotypical segmentation results.443.1. Clinical Target Volume Estimation from TRUS Images(a) (b)Figure 3.4: a) Histogram of the volumetric error and surface distance mea-sure. b) TRUS images of the marked outlier (in red color) from base (topleft) to apex (bottom right). The prostate boundary near base and apex isnot visible.3.1.4 ConclusionBy using joint independent component analysis on a spatially ordered jointdata matrix of TRUS images and their corresponding CTV labels, we haveintroduced a fast segmentation approach that can be potentially used inreal-time brachytherapy dosimetry applications. The proposed segmenta-tion approach generates a set of joint image-label basis vectors. Coefficientsobtained from partial projection of the intensity image onto these basis vec-tors is further used to combine the label basis vectors to produce a label mapfor the target image. Evaluation of this method on a clinical brachytherapydataset shows results with clinically acceptable accuracy and repeatabilitywhile it performs in less than two seconds. Accordingly, an optimized im-plementation of the algorithm with a larger dataset has the potential tohandle the requirement of real-time dosimetry for brachytherapy with intra-operative planning.453.2. Planning Target Volume and Minimum Prescribed Isodose Estimation from TRUS Images3.2 Planning Target Volume and MinimumPrescribed Isodose Estimation from TRUSImages3.2.1 IntroductionEssentially, PTV is the reference anatomical target in the preplanning ofprostate brachytherapy [5, 6, 83]. As explained in Chapter 1, the PTV isobtained by inhomogeneous expansion of the CTV according to the insti-tutional guidelines and the prostate anatomical regions. Subsequently, thenumber and pattern of seeds and needles is determined. Seed planning re-quires adherence to selected control parameters, such as the percentage ofPTV receiving the minimum Prescribed Dose value (mPD).Here, we aim to incorporate prior knowledge in a dataset of treatmentrecords including TRUS volumes, PTVs and the seed coordinates to simul-taneously predict the PTV and the isodose contour for a new patient, solelybased on the TRUS volume. The isodose contour of the prescribed dose canbe considered as one of the plan representation parameters, as it reflects themPD distribution inside and outside of the PTV. We calculate the isodosecontours based on the seed coordinates and adapt the joint IndependentComponent Analysis (jICA) [78–81], as a reconstruction based modelingapproach to capture relation between TRUS intensity volumes, PTV andthe prescribed isodose map. We propose a model with a set of joint compo-nents that are identified using jICA. The PTV and the isodose map of a newpatient TRUS volume are represented as a weighted sum of the identifiedjoint components in the derived model.A major shortcoming of the ICA-based modeling is the instability of theidentified sources, since majority of the mathematical solutions proposedfor ICA are sensitive to initial conditions. To alleviate this issue, we repeatthe ICA decomposition process multiple times. Stable components are thenchosen from the centroids of the clusters representing compactness and sep-arability of the decomposed space. In this chapter, we present a stabilityanalysis for the jICA approach. Subsequently, we follow a holdout valida-tion approach to evaluate our generated model from fusion of TRUS, PTVand the prescribed isodose. We report the result of our analysis on dataobtained from 120 patients.This section is adapted from [82]: S. Nouranian, M. Ramezani, S. S. Mahdavi, I.Spadinger, W. J. Morris, S. E. Salcudean, and P. Abolmaesumi, Data fusion for planningtarget volume and isodose prediction in prostate brachytherapy, Proc. SPIE, vol. 9415.p. 94151I-94151I-7, 2015.463.2. Planning Target Volume and Minimum Prescribed Isodose Estimation from TRUS ImagesObservation MatrixJoint 3-modality ComponentsCoeffDataset of Treatment KnowledgeUnseen VolumePTV MapIsodose MapPTV EstimatemPD EstimateStable Joint Componentsxv xl xdKCdCv ClFigure 3.5: Isodose prediction pipeline consists of a modeling (offline) pro-cess and a prediction (real-time) process. Prediction process simultaneouslyestimates PTV and isodose contours solely based on an unseen TRUS vol-ume.3.2.2 MethodThe block diagram of the proposed method is shown in Figure 3.5. A ret-rospective dataset of treatment records is used to identify joint componentsof the ICA model. The dataset consists of three information modalities;namely, the TRUS volume, PTV binary volume and the isodose binary vol-ume which is obtained by performing dosimetry calculations on the seedplan. The modeling procedure takes place in an offline mode. The gener-ated model encodes the information within the dataset into a set of jointcomponents by performing the joint component stability analysis. In theprediction phase, an unseen TRUS volume is only partially projected ontothe TRUS part of the joint components to obtain the appropriate coeffi-cients that minimize the reconstruction error. Same coefficients are adaptedto combine the other modality parts of the joint components and generatea probability map of the corresponding PTV and isodose. We apply a maskto convert these maps into the target binary PTV and isodose estimations.We denote a TRUS volume as V (p) ∈ R, the associated PTV volumeas S(p) ∈ {0, 1}, and the corresponding mPD volume as D(p) ∈ {0, 1} in a473.2. Planning Target Volume and Minimum Prescribed Isodose Estimation from TRUS Imagesdiscrete space p ∈ Ω, where Ω ⊂ R3. A unique spatial pattern u, in spaceΩ, is used to convert the original volumes into row vectors xv, xl, and xd,respectively, so that each voxel is represented with the same index in allthree vectorized modalities. We form a joint observation vector [xvxlxd] byconcatenating these data vectors, and hereafter refer to ith patient recordwith x(i).mPD Volume GenerationAccording to the guidelines used at our local cancer centre, we use theupdated dosimetry formalism proposed by Task Group #43 [84] in orderto obtain the target isodose contours from seed locations. We follow thepoint source (1-D) formalism for the radio-active seeds. The general 1-Dformalism introduces minimal error and simplifies the dose calculation byignoring the seed orientation. The maximum initial dose D˙ is expressed asD˙(r) = SK .Λ.(r0r)2.gp(r).φan(r), (3.3)where SK represents air-kerma strength of the radio-active seed (in µGym2h−1,aka U), Λ is dose-rate constant in water (in cGyh−1U−1), r0 denotes thereference distance which is 1 cm from the seed center, gp accounts for pho-ton attenuation and scatter effects and φan accounts for anisotropies due toseed construction. Total dose D deposited by a single point model source iscalculated asD(r) = D˙(r)∫ t=Tt=0e−λtdt, (3.4)where T is the time of radiation (in h). The isodose contour of the prescribeddose (usually set at 144 Gy) outlines the tissue volume that receives theprescribed mPD or greater. Maximum PTV coverage and conformity of thisdose distribution are important factors considered in devising a treatmentplan for a new patient.Model GenerationA joint three-modality observation matrix X is formed by stacking all patientrecords x(i), i = 1, 2, ..., N , where N is the size of training dataset. Wegenerate a set of K joint independent components c(1), c(2), ..., c(K) using483.2. Planning Target Volume and Minimum Prescribed Isodose Estimation from TRUS Imagesthe InfoMax approach [85], so thatC = WXˆ, where C =c(1)v c(1)l c(1)dc(2)v c(2)l c(2)d.........c(K)v c(K)l c(K)d , (3.5)where W is the demixing matrix and each component c(j) is divided intothree modality parts [c(j)v c(j)l c(j)d ], and Xˆ is the whitened observation ma-trix. The ICA algorithm is sensitive to the whitening process which involvescentering the data vectors and transforming them so that the covariance ma-trix is close to identity. The whitening process decorrelates the data vectorsprior to the decomposition step. In general, ICA works with higher orderstatistics of the observation data and assumes that the hidden componentsare statistically independent, non-Gaussian with linear mixing process.Component StabilityThe ICA decomposition algorithms do not necessarily identify stable andreliable components from an existing dataset. This uncertainty correlateswith the inherent properties of the ICA contrast function and the optimiza-tion techniques undertaken. Moreover, the observation dataset is mostlylimited in size and may induce statistical errors in the component estima-tion process. One of the approaches to overcome this problem is to repeatthe decomposition process multiple times and introduce an average estima-tion of the results. Similar to the approach proposed in [86], we repeat thejoint ICA decomposition m times in order to obtain the K components.Demixing matrices of Wi are stacked into a matrix W˙ in order to calculatethe mutual correlation coefficient fromR = W˙XˆXˆTW˙T , (3.6)where an element of matrix R represented by rij refers to the mutual cor-relation between components c(x) and c(y) in two different decompositionattempts. A hierarchical clustering on the dissimilarity coefficients 1− |rij |groups the corresponding components in different decomposition attemptsinto clusters. We introduce centroids of the clusters as the stable joint com-ponents of the model.A key parameter in our modeling approach is the number of joint com-ponents, K. Ideally, minimum reconstruction error is achieved when the493.2. Planning Target Volume and Minimum Prescribed Isodose Estimation from TRUS Imagesmaximum possible number of components is considered in the ICA algo-rithm. However, some of the generated components will reflect the system-atic noise and the variability enforced by the joint decomposition algorithm.A trade off between noise cancellation and reconstruction error minimizationis targeted by a hold-out validation approach. We randomly separate M pa-tient records for the validation dataset. The validation set is never used inthe prediction experiments. We choose the optimal number of componentsby observing the reconstruction error in three modalities on the validationdataset.PredictionFor a new unseen TRUS volume xtv, we solve a least squares problem to findthe projection coefficients a, so thatxtv ≈ aTc(1)vc(2)v...c(K)v . (3.7)According to the modeling constraint, same coefficients can be adaptedto reconstruct the other two modalities, the PTV and the mPD volume.However, a linear weighted combination of the PTV and isodose parts ofthe components produces probabilistic estimates x˜tl and x˜tl for the targetbinary volumes. We propose to generate a threshold mask from the trainingdataset used for modeling. For each 2-D axial plane of the sparse volumeout of S planes, we calculate optimal plane threshold τˆs, s = 1, 2, ..., S:τˆs = argminτ1NN∑i=1|x(i)l − b(aTi Cl, τ)|, 0 < τ < 1, (3.8)where b(vT , τ) is a function that converts the normalized scalar volume vinto a binary volume by applying the threshold value τ . For any test volume,after obtaining the probability segmentation map, we apply τˆs to plane s, inorder to estimate the target segmentation contour in that plane. Further-more, edge contours of the estimated binary volume are smoothed slice byslice to produce the final estimation of the target segmentation. Contoursare represented in a polar coordinate system assuming their geometrical cen-ter as the origin and further transformed into the Euler coordinate systemto perform B-spline curve fitting. The generated contours are smoothed toproduce the final estimation of the target segmentation.503.2. Planning Target Volume and Minimum Prescribed Isodose Estimation from TRUS Images3.2.3 ExperimentsRecords from a cohort of 120 patients who had undergone brachytherapytreatment were used in the proposed modeling approach. We only includedTRUS volumes of axially approximately equal length prostates in the cohortof this study. We enforce this constraint to perform the feasibility studywhile satisfying the approximate anatomy alignment condition. For eachTRUS volume, the PTV and seed coordinates are provided by expert radi-ation oncology practitioners. We calculate the dose distribution map andsubsequently, the mPD volume of 144 Gy based on a point source modeldescribed in Section 3.2.2. TRUS volumes consist of nine equally spaced(5 mm) 2D ultrasound images, collected from base to apex of the prostatein the transverse plane and in accordance with the VCC protocol.In all analyses, we randomly use 80% of the cases in our training datasetand evaluate the hypothesis on the remaining 20%. To avoid bias in themodeling process, we repeat the jICA decomposition with stable componentsfor 10 times.Predicted PTV and mPD volumes are evaluated against the gold-standard.We evaluate PTV with volume error (Verr) and volume difference (Vdiff ) asshape and size dissimilarity metrics [1, 54, 60]. In addition to the volumetricmeasures, we validate the predicted mPD volume based on the dosimetryparameter V100. The V100 indicates the percentage of the PTV volumereceiving 100% of the prescribed dose.3.2.4 Results and DiscussionTo determine the optimum number of components in our modeling approach,we calculate the reconstruction error for PTV and mPD volume using differ-ent numbers of components as shown in Figure 3.6. Intuitively, the numberof components in our model represents the complexity of the data vectorspace as well as the power of the model in capturing the correlation betweenmodalities. The smaller the number of components, the larger the recon-struction error observed due to the fact that the model is unable to capturethe correlation between modalities. As the number of components grows,the reconstruction error is reduced; however, there is an increase when thenumber of components is close to the total number of cases included inthe training dataset. This can be due to the stronger presence of compo-nents that reflect noise in the training dataset. Our linear decompositionapproach equivalently optimizes the weights for all components, hence isprone to miss-guidance by the these components. We choose the optimum513.2. Planning Target Volume and Minimum Prescribed Isodose Estimation from TRUS ImagesFigure 3.6: Reconstruction error Verr calculated for PTV and Isodose ontraining and validation datasets.00.050.10.150.20.251 2 3 4 5 6 7 8 9 10PTV volume error00.050.10.150.20.251 2 3 4 5 6 7 8 9 10mPD volume errorFigure 3.7: Verr statistics shown with boxplots from 10 random datasetpartitioning into train and test sets.number of 40 components as a compromise on lower reconstruction erroramong the three modalities.Figure 3.7 shows the statistics of the Verr in 10 random dataset parti-tioning. The obtained results show that PTV and mPD volume predictionerrors are not statistically different. It can be inferred that the ICA-basedmodeling approach is able to equally capture the correlation between bothPTV and isodose volumes and the original TRUS volume.The calculated grand mean and pooled standard deviation of the metricsfrom 10 trials are shown in Table 3.2. We calculated V 100 with respect tothe gold standard PTV. At the Vancouver Cancer Centre, clinical guide-lines recommend a value greater than 98% for the dosimetry parameter ofV 100. If we include the predicted PTV in V 100 calculation, we observe thatthe proposed modeling approach always meets the planning standard V 100measure, i.e. 98.46 ± 1.09% ≈ 98% compared to the gold standard V 100measured at 97.51 ± 0.56%. Paired t-test statistical significance analysisfails to show any difference in the mean value of the calculated V 100 for523.2. Planning Target Volume and Minimum Prescribed Isodose Estimation from TRUS Imagesthe proposed method and the gold-standard with p-value< 0.05. Five 2DTRUS images (from base to apex 1 cm apart) of a randomly chosen patientrecord are shown in Figure 3.8. Top row shows the treatment gold standardcompared to the bottom row as the predicted PTV and mPD volume maskoverlaid on the TRUS images.Table 3.2: Grand mean and pooled standard deviation of Verr, Vdiff andV 100 calculated for the estimated PTV and mPD volume in 10 differentrandom training/testing datasets.PTV IsodoseVerr Vdiff Verr Vdiff V10010.02± 4.5% 7.55± 14.77% 9.74± 4.23% 9.21± 14.17% 97± 3.55%3.2.5 ConclusionWe proposed a reconstruction based modeling approach to find the correla-tion between TRUS volume, PTV and the mPD volume. We employed thejICA algorithm on a set of treatment records in order to introduce joint inde-pendent components representing the correlation between these modalities.Instability of the components caused by the jICA algorithm was reviewed byclustering the components based on the correlation between jICA demixingcoefficients in different trials. We also introduced parameter tuning proce-dures for the number of components, and the binary conversion mask forthe mixed components.The proposed approach enables us to project the prior knowledge in aset of brachytherapy treatment records onto an unseen TRUS volume. Con-sequently, it can lead to reduced inter- and intra-observer variability associ-ated with the brachytherapy PTV delineation and seed planning processes;as the dose plan and expert’s knowledge can contribute to the treatmentoutcome and morbidity rate [87, 88]. The proposed modeling approach isapplicable in risk stratification and toxicity/morbidity prediction. It runsin less than a couple of seconds on a typical PC and learns the correlationpattern inside the dataset regardless of the prostate boundary and texturefeatures, therefore, is an applicable approach to capture institution/clinicianspecific knowledge. Future work will include seeds coordinates determina-tion in addition to the 150% isodose contours which are not contiguous inall planes.533.2. Planning Target Volume and Minimum Prescribed Isodose Estimation from TRUS ImagesFigure 3.8: Gold standard (top) and predicted (bottom) TRUS, PTV andmPD volume (base to apex, from left to right, 1 cm apart) of a randomlychosen test case.54Chapter 4Sparsity-based Fusion forBrachytherapy Preplanning4.1 IntroductionThere is no evidence supporting the idea that joint patterns between volu-metric information of brachytherapy such as TRUS and CTV volumes arestatistically independent. Although, this assumption has led the model-ing approach to look for joint patterns that are spatially locally distributedin the extent of the joint space. Hence, sparse pattens can potentially be anew perspective in defining the joint patterns in the fused information space.Moreover, sparsity can be seen as a generalization of the ICA-based model-ing that produces joint components that are sparse instead of statisticallyindependent. Hence, here in this chapter we investigate sparse analysis andits strength in representing the joint patterns in the fusion framework.In this chapter we aim to propose a fusion algorithm that automatesdelineation of both CTV and PTV contours. Although there are some clin-ical guidelines but definition of the CTV and the PTV is not definite andvaries from case to case. Hence, most of the previous TRUS segmenta-tion approaches are not directly applicable to CTV and PTV segmentation.We have visualized two sample actual treatment cases in Figure 4.1 to em-phasize on the difference between prostate boundary and the two targetvolumes. Figure 4.1 shows some anatomical structures in the anterior re-gion of the mid-gland that are outside the prostate boundary but insidethe CTV contours. In general, CTV is not necessarily the real prostateboundary and is defined for each patient to maximize the likelihood of re-ceiving prescribed radiation dose for the target anatomy while minimizingthe irradiation to the surrounding tissue. According to International Com-mission of Radiation Units and Measurements (ICRU Report 50) and ABS,This chapter is adapted from [89]: S. Nouranian, M. Ramezani, I. Spadinger, W. J.Morris, S. E. Salcudean, and P. Abolmaesumi, Learning-based Multi-label Segmentation ofTransrectal Ultrasound Images for Prostate Brachytherapy, IEEE Trans. Med. Imaging,vol. 35, no. 3, pp. 921-932, 2016.554.1. IntroductionBase / Base + 0.5 cm Mid-gland Apex / Apex + 0.5 cm Patient 1 Patient 2 Original: With contours: Original: With contours: PTV CTV ? Figure 4.1: Sample CTV and PTV contours used for treatment of two pa-tients. CTV is not necessarily the prostate boundary: 1) the prostate bound-ary is not visible in near base and apex views (patient 1 base views), and2) in mid-gland views some surrounding tissue is included in CTV (patient2 mid-gland and apex views). PTV-CTV margin varies within and betweenpatients (patient 1 mid-gland views); hence, PTV is anisotropic dilationof CTV per patient based on CTV, guidelines and expertise of radiationtherapist.564.1. IntroductionCTV is an anatomical-clinical concept that contains cancerous tissue andsub-clinical microscopic malignant disease, and may contain adjacent tissuessuch as extra-capsular extensions of the prostate. Although there are someguidelines for delineation of the PTV based on CTV, PTV is also deter-mined based on some patient specific clinical parameters (prostate length,size, shape and visibility) and more specifically the clinician’s expertise. Avisual comparison of the two TRUS images in Figure 4.1 shows that PTVcontours are not a predefined dilation of the CTV and size of the marginin different regions of the prostate varies from case to case. Moreover, theprostate boundary is not usually visible close to the base and the apex ofthe prostate in TRUS axial views. Furthermore, most of the previous TRUSsegmentation reports either need user interaction or require access to denseTRUS volumes. Also our approach for multi-atlas-based segmentation (seeChapter 2) shows clinically acceptable results [54, 60, 61], but none of thesemethods meet requirements for a real-time, automatic segmentation appli-cation such as speed which is highly demanded for intraoperative clinicalapplications. Moreover, none of the methods provide two target anatomi-cal labels for TRUS voxels, i.e. CTV and PTV, and therefore require twoindependent, unconstrained and sequential segmentation attempts.In this work, we aim to incorporate prior knowledge in a large datasetof treatment record observations including sparse TRUS volumes and theCTV/PTV labels to introduce a fast and automatic multi-label segmenta-tion method. We introduce our solution in the form of a joint model thatcaptures the relation between the three observations. The proposed modelis used to simultaneously estimate the CTV and the PTV labels for a newpatient, solely based on the TRUS volume. We propose a sparse representa-tion of joint observations and train a set of joint sparse dictionary elements,aka atoms. Each atom is a non-linear model of a particular joint patternbetween three volumetric information spaces. Consequently, projection of apartial observation (i.e. TRUS volume) onto the trained dictionary elements(i.e. joint atoms) estimates a set of coefficients suitable to combine otherparts of the atoms that correspond to CTV and PTV labels.To the best of our knowledge, this work is the first framework proposedfor fast and automatic prostate brachytherapy treatment target labeling ofthe TRUS images. In summary, key contributions of the proposed methodare: 1) We constrain the dictionary learning algorithm to not only satisfythe sparsity but also to force equal contributions for each volumetric in-formation space of the joint atoms. Hence, we obtain an optimum set ofcoefficients from one volumetric information space (i.e. TRUS) and apply itto the other expected volumetric information spaces (i.e. CTV and PTV).574.2. Methods2) In contrast to other statistical shape modeling approaches, this methodis not iterative (fast) and does not rely on local shape and appearance char-acteristics and incorporates all voxel values. Architecture of the proposedframework including choice of the pre-processing and the proposed thresh-olding algorithm to convert estimated probability maps of CTV and PTVinto binary labels are specifically tailored to the context of this study.We evaluate the proposed method on 590 patient records and compare es-timated contours against the manually segmented treatment target contoursthat were used by clinicians to plan the seed arrangements. A 5-fold crossvalidation shows clinically acceptable results that also satisfies speed require-ments for a real-time clinical application. We compare performance of theproposed multi-label model with single-label models generated from jointTRUS/CTV and joint TRUS/PTV labels, and show advantage of multi-label model over single-label models. We also show an added value of theproposed method: it can adapt to different clinics, clinicians or treatmentgroups, as it inherently captures the expertise and implicit guidelines usedin different centres in the form of a dataset specific sparse dictionary.4.2 MethodsFigure 4.2 illustrates the proposed pipeline for multi-label segmentation ofTRUS images. We propose a joint sparse dictionary learning approach ina one-time offline process of learning performed on a set of observationsfrom TRUS, CTV and PTV. Subsequently, the generated model is used ina real-time process of estimation. In the following sections, we introduceobservation space, modeling and its successive predictors.4.2.1 Problem Definition and Observation Space FormationWe represent patient volumetric information in a unique discrete Cartesianspace of Ω ⊂ R3. We assume all patients’ volumetric information datasetsare coregistered in this Cartesian space for further analysis. We define threevolumetric information spaces of target anatomy within the discrete spaceΩ as TRUS intensity and corresponding binary labels for the CTV andthe PTV. We assume there exist certain observable patterns for each voxelw.r.t. the other voxels within Ω in and among three volumetric informationspaces. Our goal is to learn joint patterns from a dataset of observationsby including multiple volumetric information spaces in one unique learningalgorithm. Our approach is to generate a joint sparse dictionary of ele-ments [90], i.e. atoms, so that each volumetric information space can be584.2. MethodsProjection Learning (offline) Estimation (online) Target TRUS ~ PTV Smoothing Confidence Map TRUS PTV Same Spatial Order = Smoothing Reconstruction Thresholding T Weights X CTV Dictionary Atoms T TRUS PTV Label CTV Label ~ CTV … … … … … … Figure 4.2: Illustration of the proposed system for simultaneous automaticlabeling of the TRUS volumes.represented by a weighted linear combination of these atoms. For each pa-tient i ∈ {1, 2, ..., N} in a dataset of N patient treatment records, we use aglobal spatial indexing order in the space Ω to generate three observationvectors of v(i) ∈ Rd×1, l(i)ctv ∈ {0, 1}d×1 and l(i)ptv ∈ {0, 1}d×1, representingTRUS intensities, CTV and PTV labels, respectively, where d is the num-ber of discrete coordinates in Ω. Therefore, each voxel is represented withthe same index in all vectorized triplets of volumetric information spaces.We introduce a full joint observation vector as concatenation of the threeobservations [v(i)T, l(i)Tctv , l(i)Tptv ]T .Presence of noise in observation vectors can potentially decrease perfor-mance of the training process and consequently make the atoms less likelyto capture the correlation among input volumetric information spaces. Thisis important when it comes to the noisy and speckled TRUS images. Weaim to reduce the input noise in TRUS intensities by adapting a Gaussiansmoothing kernel to convolve with 2-D axial images. The filtered images arethen used to generate smoothed TRUS observation vector of vˆ(i).CTV and PTV delineations are also affected by both intrinsic noise andobserver uncertainty. Hence, we adapt an approach [91] to provide confi-dence measures from binary segmentations. We aim to assign a confidencedegree to voxels based on their distance to the segmentation boundary. Orig-inal CTV and PTV binary volumes are converted into signed distance maps,and vectorized to produce signed distance vectors r(i). Then we apply a sig-moid function to calculate voxel’s likelihood to be inside the target label594.2. Methodsbinary volume:rˆ(i) = P (r(i)) =11 + exp(−αr(i)) , i ∈ 1, 2, ..., N. (4.1)This transformation produces a likelihood map, in which the points in theneighborhood of the target boundary are associated with lower confidencedegree than the points sitting at the boundary. Parameter α, the slope ofthe sigmoid function, determines the neighborhood size for confidence degreecalculation. This parameters is empirically set to 0.8 to represent ∼ 3 mmmargin of less than 95% confidence in label assignment.Finally, we generate the joint observation matrix Xˆ by stacking pro-cessed joint observation vectors of x(i) = [vˆ(i)T , rˆ(i)Tctv , rˆ(i)Tptv ]T . Each voxelcan be considered as a feature in the learning process; however, due to thestationarity characteristics of the observations from voxels, features in aneighborhood reflect a high correlation. We decorrelate the observation fea-tures by whitening the input observation matrix. Hence, we train the modelon the whitened observation matrix X˜ to obtain a joint sparse dictionaryrepresenting global joint patterns in Ω for the three volumetric informationspaces.4.2.2 Joint Sparse Dictionary LearningIn our joint learning approach, we aim to introduce D ∈ R3d×K as a dic-tionary of K atoms that satisfies ‖x(i) − Dc(i)‖ ≤ , where c(i) ∈ RK×1is a vector of weights obtained for joint observation vector x(i) with fewestnumber of nonzero elements.Since we introduce the joint observation vector x(i) to the learning algo-rithm, the generated dictionary D contains atoms that jointly represent theobservations in three parts, i.e. D = [DvT ,DlctvT,DlptvT]T . Hence, atomcoefficients are constrained to be equal for all three volumetric informationspaces while satisfying the sparsity and optimal reconstruction conditionsthroughout the learning process.We reformulate the training process for the whitened joint observationmatrix X˜ ∈ R3d×N and express the joint dictionary learning as X˜ ≈ DC,where C ∈ RK×N refers to the coefficient matrix [92].We use the K-singular value decomposition (K-SVD) method [93] totrain the joint model and generate the dictionary matrix D. K-SVD is ageneralization of the K-means clustering algorithm for vector quantization,604.2. Methodsand can be defined as the following optimization problem:minD,C{‖X˜−DC‖2F}subject to ∀i, ‖c(i)‖0 ≤ T0, (4.2)where ‖.‖F is the Frobenius matrix norm and T0 is the sparsity level or thenumber of nonzero coefficients. K-SVD as a minimization problem solverconsists of two steps: 1) assuming D is fixed, an orthogonal mapping pursuitalgorithm [94] is used for sparse coding, i.e. obtaining c(i) by decouplingthe optimization problem into N individual sub-problems. 2) Dictionaryatoms are updated one by one. Assuming D and C are fixed except onecolumn (joint atom) in the dictionary matrix, i.e. dj = [dvjTdlctvjTdlptvjT]T ,an optimization problem is defined with a penalty term as:‖X˜−DC‖2F = ‖Ej − djcj‖2F , (4.3)where Ej is the summation of the error for all N observations when jthatom is removed and cj is a vector of coefficients of that joint atom forall observations. Then, representation error Eˆj is calculated from thoseerror columns corresponding to observation that use atom dj (with nonzerocoefficient). Eventually, atoms and their corresponding sparse weights aresequentially updated using the singular value decomposition:Eˆj = U∆VT (4.4)Solution for dj is then defined as the first column of U and coefficient vectorcj as first column of V multiplied by ∆(1, 1).Since the K-SVD approach is sensitive to the initial conditions, i.e. theinitial dictionary, we reduce this sensitivity by conducting independent com-ponent analysis on training dataset. We adopt the InfoMax approach [85]to obtain K initial dictionary atoms for the joint sparse dictionary learningalgorithm. We let the K-SVD algorithm iterate through both steps for 30times to generate a sparse joint dictionary D in our offline learning process.4.2.3 Multi-label EstimationIn a real-time estimation process (see Figure 4.2), a new unseen TRUS imageset vt is smoothed by the same Gaussian kernel to provide an incomplete ob-servation vector of vˆt. We solve a least squares problem, i.e. the projectionstep, to find the optimal set of dictionary coefficients ct asc˜t = argminct‖vˆt −Dvct‖22, (4.5)614.3. Materialswhere ‖.‖22 denotes the Euclidean norm.We use the obtained coefficients from the projection to linearly combinethe other two parts of the atoms (related to CTV and PTV observations)and calculate P (ltctv) = Dlctvct and P (ltptv) = Dlptvct. However, P (ltctv) andP (ltptv) are probabilistic estimates of the CTV and the PTV for the test case.We propose to generate a threshold mask from the observation matrix Xˆto convert the probability estimations into a binary estimation of the CTVand the PTV. We calculate optimal plane-wise threshold τˆq, q = 1, 2, ..., Qfor CTV and PTV labels separately:τˆq = argminτ1NN∑i=1|l(i) − δ(Dlc(i), τ)|, 0 < τ < 1, (4.6)where Q is the number of axial planes acquired for the volume study in spaceΩ, and δ(y, τ) is a threshold function that converts the normalized scalarobservation vector of y into a binary vector by applying the threshold valueτ . We apply τˆq to plane q for the test case, in order to estimate the CTVand the PTV labels per plane.4.3 MaterialsA retrospective cohort of 590 patients who underwent prostate brachyther-apy at the Vancouver Cancer Centre (VCC) is used to evaluate the proposedmethod. For each patient, a complete set of observations consisting of TRUSimages, CTV and PTV contours is collected. Volume study routine in VCCinvolves TRUS image collection from the base to the apex of the prostatecaptured every 5 mm in the axial direction using a side-firing endorectalprobe, hence the resulting image set is a sparse representation of the volumeof the prostate and the adjacent tissue. According to the size of the prostate,total 7-14 TRUS images of size 415×490 pixels (physical spacing of 0.1557mm by 0.1560 mm) were collected per patient. The institutional clinical vol-ume study protocol obliges the prostate to be maintained in the middle ofthe axial plane, making the anatomy appear symmetric or nearly symmetricin all ultrasound images. These series of 2-D images are then contouredsemi-automatically [1] and further modified to produce CTV contours. Fol-lowing the clinical guidelines and radiation oncologist expertise, the PTVcontours are drawn. We convert both contour sets into binary volumes forthe purpose of our study. According to the aforementioned workflow, TRUSvolume, CTV and PTV are all confined axially within the base to the apex624.4. Experiments and Resultsof the prostate for each patient. Therefore, target anatomy is implicitly co-registered among all patients. No re-sampling of the images was necessary,since, all TRUS volumes were acquired with a consistent depth parameterand the same transducer type; hence, all pixel observations reside on theunified coordinate system of Ω.We evaluate the performance of the multi-label estimator and compare itwith contours of the actual delivered treatment plans, by calculating volumeshape and size similarity coefficients; namely, Verr and Vdiff in differentanatomical sectors of the prostate [24]. We also measure distance betweenactual and estimated CTV and PTV surfaces measures in terms of HausdorffDistance (HD) and Mean Surface Distance (MSD).4.4 Experiments and ResultsA 5-fold cross validation approach was implemented to measure accuracyof the multi-label estimation. We performed some experiments to tune thehyper-parameters of the proposed model. Subsequently, we validated themodel on all the 590 cases to review its performance in general. Addition-ally, we divided the cohort into physician specific sub-groups in order toinvestigate adaptability of the method. In all experiments we implementeda 5-fold cross validation approach after hyper-parameter tuning.4.4.1 Model Parameter TuningWe randomly partitioned the cohort into 80% training and 20% validationgroups, to tune the hyper-parameters of the proposed model, namely, num-ber of atoms and level of sparsity. Partitioning for parameter tuning isperformed once to obtain an estimation of the complexity of the joint infor-mation in the dataset. Although tuning can be based on multiple observa-tions from repeated random partitioning, we did not find performance of theproposed algorithm be very sensitive to these parameters. We characterizedperformance of the model by measuring volume errors for estimated CTVsV ctverr on the validation dataset. An exhaustive grid search is carried outto determine complexity of the patterns in terms of the minimum numberof atoms and sparsity level to estimate CTV labels from TRUS volumes.Figure 4.3 depicts mean and standard deviation of CTV Verr calculated fordifferent values of sparsity level and dictionary size. Cold region (shownin dark blue) in both heat maps represents combination of the two hyper-parameters that produces results with both lower mean and lower standarddeviation of the error. It can be observed that increasing dictionary size634.4. Experiments and ResultsVerr MeanSparsity LevelDictionary Size 20 40 60 801020304050 1012141618Verr Standard DeviationSparsity LevelDictionary Size 20 40 60 80102030405044.555.566.57(a)Verr MeanSparsity LevelDictionary Size 20 40 60 801020304050 1012141618Verr Standard DeviationSparsity LevelDictionary Size 20 40 60 80102030405044.555.566.57(b)Figure 4.3: Hyper-parameter tuning by grid search on sparsity level anddictionary size; a) the mean, and b) the standard deviation of CTV Verrcalculated for validation group, i.e. random partitioning of the dataset into80% training and 20% validation. Dashed line shows the optimum parame-ters used in this study.and sparsity levels includes larger number of atoms that represent noise inthe system. Moreover, very small sparsity level fails to learn patterns in thedataset. As a trade-off between estimation performance and fewer numberof atoms and coefficients, we chose a dictionary size of 60 atoms with spar-sity level of 11, as we observed lower validation error variation in estimatingthe CTV.4.4.2 Multi-label EstimationFigure 4.4 illustrates some of the patterns captured by the proposed jointsparse dictionary learning. It is observed that atoms jointly represent spa-tially distributed patterns as is highlighted in Figure 4.4. Each atom shownin this figure represents observation variations in different regions, e.g. base-anterior, posterior, lateral and apex-anterior. This figure also shows howjoint patterns are captured, as the PTV patterns are highly correlated withCTV patterns per each atom.General statistics of the estimation performance characterized by volu-metric measures in different anatomical sectors of the prostate are shown inTable 4.1 and Figure 4.5. We also measured MSD and HD as 0.95 ± 0.41mm and 5.33± 1.47 mm for CTV, and 1.19± 0.48 mm and 5.48± 1.51 mmfor PTV.It can be observed that a lower estimation error is obtained in mid-644.4. Experiments and ResultsFigure 4.4: Learned sparse dictionary atoms. Each joint atom representsregion-specific modulation in base-anterior, posterior, lateral and apex-anterior regions for all three volumetric information spaces. Highlightedarea shows coordinates with more than 50% modulation value within thesame volumetric information space.Table 4.1: Mean and standard deviation of Verr calculated for nine anatom-ical sectors of the prostate over CTV and PTV.Verr% Base Mid Apex TotalCTV PTV CTV PTV CTV PTV CTV PTVAnt. 11.33±8.53 10.79±7.63 8.75±6.91 8.38±6.36 13.81±9.46 12.69±8.839.92±3.51 8.84±3.13Lat. 8.42±4.19 7.40±3.32 7.41±3.53 6.49±3.15 11.35±5.48 9.59±4.76Post. 8.87±6.79 7.91±6.06 7.44±6.31 7.06±5.97 12.15±9.05 11.21±8.41posterior and mid-lateral sectors of the prostate. However, in anterior sec-tors and apex region a larger deviation from the actual PTV is measured.It is reasonable to consider that the prostate boundary is poorly visible an-teriorly and particularly in close to apex and base TRUS images; hence,654.4. Experiments and ResultsPosterior Lateral Anterior Posterior Lateral Anterior Posterior Lateral Anterior01020304050607080Ve rr (%) Base Mid ApexCTVPTVFigure 4.5: Label estimation errors in terms of Verr in different anatomicalsectors of the prostate. Circles represent mean and the horizontal line themedian of the calculated error in each sector.a greater observer variability is expected. Figure 4.5 also shows that lessnumber of outliers are seen in lateral sectors. In general, we observe lowervolumetric errors for the estimated PTV in all sectors as well as the entirevolume. Two typical multi-label estimation results are shown in Figure 4.6per slice from base to apex of the prostate.4.4.3 Dataset Adaptability AnalysisWe extract three sub-groups from the cohort of 590 patient records basedon physician in charge of the treatment planning. We selected the threephysicians with the largest number of assigned cases, hereafter referred toas physician A, B and C. To review dataset adaptability of the proposedmulti-label segmentation approach, we train the model for groups A, B andC independently and perform 5-fold cross validation within each group (in-ternal validation). Subsequently, we use the trained model obtained for eachfold to validate on the remainder of the entire 590 patient cohort (externalvalidation).Hyper-parameter tuning for group A (174 cases), B (156 cases) and C(135 cases) is performed as explained in Sec. 4.4.1. General statistics of thevolumetric error measure Verr for the three groups are shown in Figure 4.7.664.4. Experiments and ResultsTable 4.2: Comparison to two other brachytherapy CTV segmentation meth-ods for sparse TRUS volumes. “?” and “•” indicate a statistically significantdifference between mean and standard deviation of the error measured bythe proposed multi-label approach and each of the reported methods usingpaired t-test and Levene’s test, respectively.Volumetric Errors (%) Surface Errors (mm)Speed (sec) AutomaticVerr Vdiff MSD HDMahdavi et al. [1] 8.81±5.60 ? • 3.56±14.66 ? 1.11±0.81 ? • 5.57±2.28 • ∼ 120− 240 NoNouranian et al. [54] 10.25±3.74 -2.37±13.73 1.33±0.60 ? • 5.48±1.71 • ∼ 180 YesProposed Method 9.95±3.53 -3.25±11.42 0.98±0.39 5.4±1.38 ∼ 0.3 Yes4.4.4 Comparison of the CTV Labeling with OtherMethodsThe proposed multi-label segmentation approach is not directly a prostatesegmentation solution. Since to the best of our knowledge, there are no priorwork on automatic PTV estimation from TRUS, we compare the perfor-mance of the CTV estimations we obtain with two of the recently publishedworks: 1) a semi-automatic approach proposed by Mahdavi et al. [1] whichhas been used at VCC for more than five years; and 2) our previous workon CTV segmentation using label fusion from multi-atlas approach [54].Current planning procedure at VCC is based on a semi-automatic ap-proach that requires manually locating some anatomical landmarks in themid-gland [1]. Contours obtained from this state-of-the-art method are fur-ther modified by clinicians to provide CTV. Then a dilated PTV volume isgenerated and manually corrected based on clinician’s expertise and knowl-edge to make sure clinical guidelines are met.Table 4.2 shows statistics of the estimated CTV error from 1) the semi-automatic algorithm, 2) the automatic multi-atlas approach, and 3) theproposed multi-label algorithm as well as the execution time. We usedthe same dataset of 280 cases that previous works are evaluated on. Weobtained 35 atoms and sparsity level of 10 by optimizing hyper-parametersof our model on this dataset, and evaluated its performance using 5-foldcross validation as explained in Sec. 4.4.A statistical significance analysis has been performed on all reported er-ror measures between the proposed multi-label approach evaluated on thedataset of 280 cases, and each of the two CTV segmentation methods usingpaired t-test and Levene’s test. As for paired t-test, we define the null hy-pothesis as equivalence between mean (t-test) and variance (Levene’s test) ofobservations from the same evaluation measure for the proposed method and674.4. Experiments and Resultseach of the other two methods. Results as indicated by “?” and “•” in Ta-ble 4.2 show statistically significant improvement for MSD when comparedto the other two methods (p-value< 0.01 for both tests). Hausdorff distanceis also significantly improved compared to the semi-automatic algorithm, byanalyzing variance of the distributions using Levene’s test. However, thesemi-automatic algorithm outperforms in terms of Verr when compared tothe proposed method. Although, according to the reported observer variabil-ity of CTV contouring, they all satisfy the clinical acceptance requirements.It is also observed that the proposed method is at least two orders of magni-tude faster than other compared methods for CTV delineation. Figure 4.8compares the proposed method with the other two CTV segmentation ap-proaches [1, 54] in terms of percentage of cases estimated with differentthreshold values for Verr. The overall observer variability is shown by avertical dashed line as the summation of reported inter- and intra-observervariabilities (see[95, 96]). The proposed multi-label method can generate ahigher percentage of successfully labeled cases w.r.t. to the overall observervariability threshold when compared to the previous work [54], but still haslower rate when compared to the semi-automatic standard-of-the-care algo-rithm proposed by Mahdavi et al. [1].4.4.5 Multi-label Model Comparison with Single-labelModelJoint dictionary atoms of the proposed method represent relation amongthree volumetric information spaces of the treatment record, i.e. TRUS,CTV and PTV. Learning process is constrained to produce equal weightsfor all three information spaces of the joint atoms. We further review sig-nificance of the relation between TRUS/CTV and TRUS/PTV by trainingmodels that generate single label from TRUS images to verify performance ofthe proposed multi-label segmentation algorithm. Three single-label jointdictionaries are trained for TRUS/CTV, TRUS/PTV and CTV/PTV tomeasure performance of the model in capturing the relation between twovolumetric components. We obtained optimum hyper-parameters of eachsingle-label model independently by performing exhaustive grid search. Weused 55 and 60 atoms for TRUS/CTV and TRUS/PTV models respectively,while sparsity level remained in the same range of 10. Performance of eachsingle-label model is compared with the proposed multi-label approach us-ing 5-fold cross validation. Figure 4.9 shows the difference between twosingle-label models and the unified multi-label model in terms of CTV andPTV Verr. Both single-label models are trained using the proposed joint684.4. Experiments and Resultssparse dictionary learning algorithm. Results show that joint patterns be-tween TRUS and CTV are better captured when compared to joint patternsbetween TRUS and PTV. This is expected as there is more variability indetermining PTV rather than CTV since PTV is determined based on CTVand the original TRUS image as well as some clinical guidelines and clini-cian’s expertise. However, if we constrain the three volumetric informationin a unified multi-label model, we achieve a higher accuracy for PTV at theprice of slightly increased error for CTV. Gray area in Figure 4.9 right, rep-resents the multi-label CTV and PTV estimation error per case. It is sortedto compare how the single-label model performs for the same case as themulti-label model. It is observed that CTV estimation is more persistent.It can be concluded that the dictionary training constrained by PTV hasminimum impact on the estimation of CTV, while in contrast, PTV estima-tion benefits from being constrained by CTV components of the observationmatrix. Our experiment on measuring performance of a two-part model tocapture correlation between CTV and PTV shows 16.28 ± 2.39% Verr inestimation of PTV from CTV. This shows a better performance when com-pared to the single-label model (TRUS/PTV) and still significantly highererror (p-value< 0.01) when compared to the multi-label approach.4.4.6 Execution TimeThe proposed method includes two independent groups of computations, oneto train the model and obtain the sparse dictionary atoms, and the otherto estimate segmentation for an unseen TRUS volume. Assuming that theinput TRUS volume is confined axially within base and apex of the prostateand is in adherence with the clinical protocols, the mathematical compu-tation cost of the estimator is trivial. The proposed method has been im-plemented in MATLABR©environment (Mathworks, MA, USA) using C++and MATLAB scripts. The computational burden of the proposed estimatoris on solving the least squares problem for TRUS projection onto the jointatom space, which is facilitated by taking advantage of efficient implementa-tion of the QR decomposition method in MATLAB environment. Therefore,we obtained very fast segmentation results from our CPU implementationon a standard PC (Intel Core i7, 2.93GHz, 8GB RAM) with execution timeless than 300 ms including smoothing pre-processing and thresholding post-processing steps. This characteristic of the proposed method makes it anappropriate choice for some real-time segmentation applications such as in-traoperative preplanning.694.5. Discussion and Conclusion4.5 Discussion and ConclusionIn this work, a fast, automatic multi-label simultaneous segmentation methodis introduced that is capable of incorporating previous treatment knowledgein the form of a joint linear model. We proposed a joint multiple volumetricinformation modeling approach that encodes relationships among multiplevolumetric information spaces, i.e. TRUS, CTV and PTV, in a set of sparsedictionary atoms. The proposed method takes advantage of stationarity ofthe observations from the anatomy jointly and in multiple volumetric infor-mation spaces. The data driven approach looks for certain joint patternsconstrained by the same projection coefficients in the form of a set of jointsources.Given the TRUS volumes are constrained with the clinical protocols, i.e.Vancouver Cancer Centre guidelines, this approach can automatically andsimultaneously provide CTV and PTV contours solely from an unseen TRUSvolume. The proposed method is resilient to ultrasound artifacts as it incor-porates voxel information w.r.t. all other voxels within the same observationcase and amongst all other observation cases. This automatic approach canpotentially be used in brachytherapy intraoperative preplanning, where theseed plans are constructed prior to the implantation procedure. This datadriven multi-label segmentation approach provides both CTV and PTV con-tours simultaneously with maximum fidelity to the observed previous targetanatomy labels, hence can produce contours that require minimum furthermanipulation in a time costly scenario such as intraoperative preplanning.4.5.1 Clinical SignificanceThe observer variability involved in locating the base and the apex axialviews of TRUS volumes as well as contouring CTV and PTV, reflects itself inthe overall labeling accuracy. It is difficult to dissociate sources of variabilityfrom other factors such as inter- and intra-observer variability of labeling,since the true prostate shape, and the base and apex planes are unknown. Inthe lack of a true gold-standard, we argue that a clinically acceptable methodshould achieve accuracy within the overall observer variability determinedfrom planning. Such analysis has been performed at VCC with severalclinicians planning for several patients [1], and inter- and intra-observervariability of manual contouring of the CTV has been reported to be in theorder of 4.65± 0.77% and 5.95± 1.59% for Verr, respectively [1]. Hence, theoverall observer variability can be estimated as the summation of these twovariabilities (see [95, 96]). We argue that this range of the overall observer704.5. Discussion and Conclusionvariability (approximately 10%) is clinically acceptable, since a statisticalanalysis of 10 years clinical outcome for VCC has shown excellent outcomeand low rate of cancer recurrence [5]. A greater variability is expected fordelineation of the PTV, since more user interaction is involved on top ofCTV labeling to determine the PTV. This potentially significant variabilityis mainly due to the low signal-to-noise ratio of the target anatomy in thesparse TRUS images, especially close to the base and the apex.4.5.2 Computational CostThis learning-based segmentation approach provides a set of joint volumetricinformation dictionary elements and hence reduces the dimension of obser-vation variations to the order of number of sparse coefficients associatedwith the atoms. Hence, the computational complexity of the multi-labelsegmentation is reduced dramatically to the order of a least squares prob-lem solution. Hence, this fast simultaneous segmentation can potentiallyfacilitate some clinical applications such as intraoperative preplanning andtargeted biopsy.4.5.3 Dataset AdaptabilityAn added advantage of the proposed data-driven multi-label segmentationapproach is the adaptability of the model to different datasets. Our experi-ments on three physician-based sub-groups within the studied cohort of pa-tients shows a high fidelity to the existing knowledge in the training dataset,as the trained model, i.e. joint sparse dictionary, performs superiorly onprediction of both segmentations for the same physician when compared tothe other physicians. Although in our study all physicians use the sameclinical protocols within VCC, our results show that the proposed methodcaptures clinician’s planning style in addition to the enforced guidelines inan institution. Hence, this method can be used in different institutions withmaximum adaptability to their current segmentation routine and guidelines.However, preplanning protocols used in this study are an evolved and ex-tended version of Seattle protocol [7] for prostate brachytherapy which iswidely used in North America since 2009. Hence, extent of the volumetricinformation spaces as well as the orientation of the prostate in axial viewsare common beyond the collected dataset in this study. These guidelinesare recommended to facilitate easy registration during patient setup; hencereduce deviation of the delivered treatments from preplans.Our proposed joint intensity and labels approach has the following unique714.5. Discussion and Conclusionkey features: 1) In contrast to the most statistical models that incorporateshape priors, the proposed method is not iterative and is super fast (com-putation time in the order of milliseconds); 2) The information encoded inthe dictionary is constrained by sparsity which provides a unique mappingthat captures spatially distributed joint patterns instead of wholistic modesof variations. Since the variability of CTV and PTV are larger in the baseand the apex in comparison to the mid-gland, we expect that spatially dis-tributed joint patterns generated by our approach can capture such localvariations with sufficient accuracy. 3) Prior intensity-shape-based models(such as [35, 36]), have only focused on prostate segmentation, hence theirapproach is not directly applicable to delineate PTV that is an anisotropicpatient-specific dilation of the CTV.4.5.4 LimitationsIn the present work, we assume volumetric information spaces of all patientsare pre-aligned according to the clinical protocols at VCC. Ideally, accuratealignment of the volumetric spaces can lead to a more accurate model (dic-tionary) or a lower number of joint components to represent the variationsamong treatment cases. In the context of our study, there are two mainlimitations that would affect the registration process: 1) Prostate shapesand sizes vary considerably; 2) TRUS image quality is low, especially in thebase and the apex of the prostate. Hence, performing rigid registration asa pre-processing requirement would introduce a potential failure point inthe overall automatic labeling workflow. Given that our results on data ob-tained at VCC are clinically acceptable without such pre-processing step [5],we conclude that enforcing a unified clinical workflow for data acquisitionwithin a center can provide TRUS data that are aligned with accuraciessufficient for our method, considering all observer variabilities.In summary, we proposed a multi-label segmentation approach evalu-ated on a retrospective dataset of treatment records at VCC. Results of theproposed models are compared against the actual CTV and PTV, and wehave shown clinically promising results that can reduce the variability asso-ciated with segmentation of the two anatomical targets. Our previous CTVsegmentation work was not directly applicable to PTV labeling since thereis low correlation between TRUS intensity and PTV labels without incor-porating CTV labels [54]. A comparison of our proposed approach to otherCTV segmentation methods show significant improvement in clinical label-ing error measures such as mean surface distance and Vdiff , and speed-upof computation time in two orders of magnitude. This method can poten-724.5. Discussion and Conclusiontially aid the clinicians by providing simultaneous estimation of CTV andPTV instantly, solely based on the collected TRUS images. The proposedsystem can be used for training, decision support or some real-time clini-cal applications such as intraoperative preplanning and targeted biopsy. Itcan be adapted to different institutions, physicians or generally the trainingdataset, and is potentially capable of predicting multiple contours that arebiased towards the implicit knowledge and expertise hidden in a datasetwhich may not be easily translated into clinical guidelines.734.5. Discussion and ConclusionFigure 4.6: Two typical results of multi-label estimation from TRUS volume.Left column in each case shows the actual and estimated (dashed line) CTVand right column the actual and estimated (dashed line) PTV from base toapex (top to bottom).744.5. Discussion and ConclusionPhysician A Physician B Physician C00.020.040.060.080.10.120.140.160.180.2CT V ( Ve rr%) Internal ValidationExternal Validation(a)Physician A Physician B Physician C00.020.040.060.080.10.120.140.160.180.2PT V ( Ve rr%) Internal ValidationExternal Validation(b)Figure 4.7: Learning-based multi-label segmentation adaptability analysis.For three different groups a model is trained and evaluated internally andexternally on the remainder of the cases.754.5. Discussion and Conclusion0 5 10 15 20 25020406080100Verr(%)T es t Ca se s ( %)CTV Estimation Proposed MethodMahdavi et al.Nouranian et al.Figure 4.8: Cumulative histogram of the CTV estimation Verr comparing theproposed multi-label approach with the two CTV segmentation algorithms:1) A semi-automatic standard of care algorithm [1]; and 2) an automaticmulti-atlas label fusion algorithm [54]. Red vertical dashed line representsoverall observer variability as a threshold to compare performance of themethods.764.5. Discussion and ConclusionCTV PTV0102030405060Verr (%) Multi-label ModelSingle-label ModelFigure 4.9: Comparison of single-label models generated for TRUS/CTVand TRUS/PTV estimation with the multi-label approach using the jointsparse dictionary learning algorithm. CTV single-label model performs bet-ter while PTV single-label model fails to capture TRUS/Label correlation.Label estimation error per case (right) shows that constraining three volu-metric information all together in the joint dictionary learning model hasminimum impact on CTV estimation (top) but a significant impact on im-proving PTV estimation when compared with single-label models.77Chapter 5Automatic Seed PlanEstimation5.1 IntroductionIn LDR prostate brachytherapy, seed arrangements can be planned eitherprior to, right before or during implantation procedure. According to theBCCA prostate brachytherapy program, preplanning is a treatment ap-proach, which requires pre-operative planning of seed arrangements, usu-ally days before the implantation procedure. Seed planning takes place inan anisotropically dilated volume around the prostate boundary, known asPlanning Target Volume (PTV). Standard-of-care for prostate brachyther-apy preplanning and delivery is to use a standard brachytherapy grid tem-plate. This template is used for intra-operative delivery of the treatment,where needles are inserted according to the preplan, and confines needleplacements to the grid points. Although clinical guidelines may vary some-what between institutions, the preplanning treatment approach has shownexcellent outcomes in the past decade [5].Currently the standard-of-care in the majority of clinics is to create theplan manually, which is subject to a certain degree of observer variability.Treatment variability also stems from large set of possible solutions that allmeet the clinical guidelines. However, a trained expert selects the recom-mended plan from this large set based on his/her expertise. Preplanningautomation according to clinical guidelines and based on previous successfulplans, can mitigate this problem.In this work, we propose a new framework to automate the prostatebrachytherapy preplanning. It consists of two components: 1) a plan esti-mator that represents sparse joint patterns between PTV and seed plans,and 2) a stochastic optimizer that encodes clinical standards to fine-tuneThis chapter is adapted from [97]: S. Nouranian, M. Ramezani, I. Spadinger, W.J. Morris, S. E. Salcudean, and P. Abolmaesumi, Automatic Prostate BrachytherapyPreplanning Using Joint Sparse Analysis, in Medical Image Computing and Computer-Assisted Intervention: MICCAI, 2015, vol. 9350, pp. 415-423.785.2. MethodsPlan Estimator Optimizer jSDL Learning (Offline) Estimation (online) Clinical Guidelines Plan Estimation Test PTV PTVs Plans Figure 5.1: Block diagram of the proposed framework. Initial estimation ofa plan is optimized, based on a novel cost function, to ensure that clinicalguidelines are met.the outcome of the plan estimator. We show that the estimator aids theoptimizer to converge to a solution that meets the clinical standards. Webuild the plan estimator by employing a sparse dictionary learning algo-rithm which encodes the joint patterns between PTV and seed plans in alarge retrospective dataset, by a set of sparse joint dictionary elements. Theestimated plan is used as an initial seed configuration for the optimizer.The stochastic optimizer is designed based on simulated annealing, whileenforcing clinical requirements through a novel cost function that capturesthose requirements. The proposed pipeline is evaluated on a dataset of 590patients who underwent brachytherapy treatment at the VCC. We demon-strate an accuracy of 86%, defined on the basis of the clinical standards andprevious treatment plans.5.2 MethodsFigure 5.1 shows the proposed framework which consists of two parts: 1) aone-time offline process of generating a model, based on the joint informationbetween PTV and seed configurations in retrospective patient data. We usea joint sparse dictionary learning approach for training the model; and,2) an online process of plan estimation and optimization. The optimizerrearranges the seed configuration in the neighborhood of the initial planestimation.5.2.1 Joint Sparse Dictionary LearningWe represent all patients’ volumetric information, consisting of the PTVand seed coordinates, in a unique discrete Cartesian space of Ω ⊂ R3, and795.2. Methodsa brachytherapy grid template space of Ψ ⊂ Ω with equal spatial spacing(e.g. 5 mm). We assume that this information is coregistered in the tem-plate space for further analysis. We define two information modalities asobservation of two sets of labels assigned to each point in the discrete spaceof Ψ, i.e. labels w.r.t. the PTV volume and labels w.r.t. the seed distri-bution in Ψ. Specifically, for each patient i ∈ {1, 2, ..., N} in a dataset ofN patient treatment records, we use a unique spatial pattern in the spaceΨ to generate l(i) ∈ Rd×1 and s(i) ∈ {0, 1}d×1, corresponding to PTV la-bels and seed placement labels in brachytherapy grid space, where d is thenumber of discrete coordinates in Ψ. We introduce joint observation vectorx(i) = [l(i)T, s(i)T]T to train a dictionary of sparse joint patterns.Our main objective is to introduce D ∈ R2d×K as a dictionary of Katoms that satisfies x(i) = Dc(i), where c(i) ∈ RK×1 is a vector of coefficientsobtained for observation vector i ∈ 1, 2, ..., N .We reformulate the dictionary training process by stacking the joint ob-servation vectors x(i) into a matrix of joint observations X ∈ R2d×N ex-pressing the joint dictionary learning as X = DC, where C ∈ RK×N refersto the coefficient matrix, and D contains atoms that jointly represent theobservations in two parts, i.e. D = [DlT,DsT ]T . Hence, atom coefficientsare constrained to be equal for both modalities while satisfying the sparsityand efficient reconstruction conditions throughout the learning process. Weuse K-singular value decomposition (K-SVD) method [93] to train the modeland generate D. K-SVD is a generalization of the K-means clustering al-gorithm for vector quantization. This is an iterative approach that consistsof two steps, 1) updating the atoms and their weights sequentially for thecolumns of D using singular value decomposition, and 2) to approximatesparsity using orthogonal mapping pursuit [94].5.2.2 Seed Plan EstimationWe perform joint sparse dictionary learning on observation vectors whichare defined on the brachytherapy grid space, i.e. Ψ ⊂ Ω. Coordinatesof the planning grid are obtained individually for each patient accordingto the clinical guidelines. Subsequently a randomly chosen case is used asreference to align all other planning grids. After transformation of PTVsand seed plans in space Ψ, we obtain aligned observation vectors of lˆ andsˆ for PTV and seed plan, respectively. We generate the observation matrixxˆ(i) = [ˆl(i), sˆ(i)] from the actual plan and the PTV to generate the joint sparsedictionary of D, where Xˆ = DC and D = [DlT,DsT ]. The dictionary Dforms the model that captures the joint relationship between PTV labels805.2. MethodsProjection Learning (offline) Estimation (online) Template Alignment Reconstruction Thresholding = X Weights Brachy Template Template Space Transformation Volume/Seeds/Needle Count Estimator Initial Estimation of the Plan Template Space Transformation Seed/Template Alignment Plan PTV Test PTV Dictionary Atoms PTV Label Seed Label T T Figure 5.2: Initial plan estimation using the joint sparse dictionary learningalgorithm.and the corresponding seed arrangement.An online estimation process includes an aligned PTV of a new case toform the observation vector lt ∈ {0, 1}d×1. Partial projection coefficientsare obtained by solving a least squares problem c˜t = argminct ‖ˆlt −Dlct‖22,where ‖.‖22 denotes the Euclidean norm. Subsequently, a probability map ofseed arrangement is calculated by linear weighted combination of the jointatoms, P (st) = Dsct.We propose to use a dynamic threshold value to convert the recon-structed P (st) into an initial estimation of the seed arrangement, s˜t. Wetake advantage of the high correlation between number of seeds and the sizeof the PTV, and generate a polynomial regression model from the trainingdataset. This model is then used to estimate the number of required seeds(according to a confidence interval), for a test case based on its PTV size.5.2.3 Plan OptimizationThe estimated seed arrangement from learning-based algorithm is not con-strained by planning quality measures, and only captures the arrangementpatterns w.r.t. the PTV coordinates; hence does not necessarily meet theclinical guidelines. We propose an optimization process to rearrange esti-mated seed configurations towards maximizing adaptation to the clinicalstandards. Table 5.1 lists some of the clinical standards that we encode into815.2. MethodsTable 5.1: Clinical guidelines at our local cancer centre:Condition GoalSeed distribution w.r.t the grid symmetricSeed distribution w.r.t the grid follows the mask150% isodose contiguousSeed spacing > 7 mmV100 > 98%V150 between 50% and 60%V200 < 21%Number of adjacent needles (row) < 3Number of adjacent needles (column) < 4the optimization algorithm. Given obtaining optimal seed arrangement is aninverse problem solution, we develop a simulated annealing (SA) algorithmthat the large search space is confined in the neighborhood of the estimatedseed arrangement s˜t.We follow the formalism proposed by Task Group #43 [84] to calculatethe dose distribution map in the discrete space of Ω. We propose an opti-mization algorithm based on SA architecture as a stochastic optimizer for aPTV volume, lptv ⊂ Ω and the corresponding estimated seed arrangements ⊂ Ψ, which we hereafter refer to as a 2-tuple y =< lptv, s >. We defineplan quality indicators as:Vp(y) : % of the PTV volume that is reached by at least p% or the mPD.Dp(y) : The minimum dose that reaches p% of the PTV volume.Cneedle(s) : Total needles required to deliver the plan.Architecture of the optimization algorithm is illustrated in Figure 5.3.We introduce a cost function that is optimized through an annealing processfollowed by a guideline-based correction algorithm. Plan is iteratively up-dated if acceptance condition is met within that iteration. A global parame-ter, aka temperature, asymptotically approaches zero as iterations progress.We use the temperature parameter in both annealing and correction pro-cesses. We define the objective function asJ(y) = αJV100(y) + βJV150(y) + γJV200(y) + δJneedles(y) + Γ(s), (5.1)825.2. MethodsCost Calculation Cool Down Anneal Seeds Apply Guidelines Acceptance Condition Update Plan Estimated Plan Optimizer Figure 5.3: Architecture of proposed optimization algorithm based on sim-ulated annealing.where α, β, γ and δ are relative importance weights assigned to differentterms of the objective function. Each term represents clinical quality in-dicators. Function Γ(y) is defined to account for contiguity of the 150%isodose by counting the non-contiguous planes and penalizing the cost func-tion. JV100 , JV150 and JV200 are introduced to calculate deviation of the dosemap from the boundary conditions:JV100(y) = 1− V100(y), (5.2)JV150(y) = u(V150(y)− 0.6)(V150(y)− 0.6) + u(0.5−V150(y))(0.5−V150(y)),(5.3)JV200(y) = u(V200(y)− 0.21)(V200(y)− 0.21), (5.4)Jneedles(y) = u(N (s)−R(y))(N (s)−R(y))R(y)+ ζ(s), (5.5)where u(x) is a Heaviside step function that returns zero for x < 0 andone for x >= 0. N (s) represents number of needles required for the seedconfiguration s, and R(y) is an estimation of the total needles expected forthe seed configuration in y. We use a high order polynomial regression modelto provide priors for needle count and PTV size relationship. ζ(s) penalizesthe cost function by returning 1, if there are more than two adjacent needlesin a row or three in a column.The SA algorithm iteratively anneals the seed arrangement s, as is ex-plained in Algo. 1. Number of seeds k and range of movement r are linearlymodified in each iteration based on the temperature cool down function. In835.2. Methodseach iteration, seeds are ranked by assigning a movement probability valuePM , which is calculated by alternating between observations of two randomvariables: 1) ratio of locally non-overlapping volume between isodose volumeand PTV, and 2) ratio of seed density over total number of seeds. Calcu-lation of PM is obtained in a neighborhood of size L which is also a linearfunction of the SA temperature. We propose the alternating approach forcalculating PM to push the annealing towards filling the uncovered areas,while making the seed arrangement more evenly distributed inside the PTV.Eventually, top k number of candidate seeds are chosen based on PM .Algorithm 1 Seeds annealing1: procedure marking candidate seeds for annealing(per iteration)2: for all seeds do3: movePrb← calculate seed move probability PM in neighborhood of size L4: candidateSeeds← top k number of seeds with highest PM5:6: procedure seed annealing7: for a randomly chosen candidate seed in candidateSeeds do8: availableSpots← find all spots on Ψ in a neighborhood of r9: for all spots do10: if spot is a forbidden spot on the grid then11: availableSpots← availableSpots - spot12:13: if availableSpots is empty then14: continue the loop.15: else16: seed← update with a randomly chosen spot in availableSpots17:18: call Seeds arrangement correction procedure Algo. 2Annealing procedure is followed by a seed arrangement correction proce-dure that is explained in Algo. 2. Seed locations, rearranged in the annealingprocess, are marked in accordance to their adherence to the clinical guide-lines, including seed spacing, brachytherapy grid mask and smaller numberof needles. Those seeds that do not satisfy the clinical guidelines are movedin a neighborhood of size r until all seeds comply with the clinical guide-lines. In our proposed SA-based optimization algorithm an exponentiallydropping temperature function with re-annealing scheme is adapted, thatgradually decreases the seed displacement range r as well as total number845.2. Methodspre-base post-apex Learning-based Initialization Arrangement Correction Arrangement Correction Annealing Figure 5.4: Seed arrangement correction and annealing process for one iter-ation of the optimization algorithm. Arrows point to a location where seedarrangement needs to be corrected based on clinical guidelines.Corrected Arrangement Local Density Map Local Overlap Map Rearranged Seeds Figure 5.5: Seeds are rearranged (annealed) based on probability valuescalculated by alternating between local seed density maps and local overlapmap between 100% isodose volume and PTV. Local neighborhood size isobtained from temperature parameter of the simulated annealing.855.2. Methodsof candidate seeds to be moved k.Figure 5.4 illustrates the annealing procedure in one iteration of theoptimization algorithm. Initial seed configuration is corrected accordingto the clinical guidelines. Then seeds are moved in a neighborhood basedon probability maps illustrated in Figure 5.5. Probability maps alternatebetween local seed density and local overlap between PTV and 100% isodosevolume, and try to fill cold regions while avoid congestion of the seeds inone region.Algorithm 2 Seeds arrangement correction1: procedure marking seeds need correction2: for all seeds do3: if seed is on a forbidden spot on the grid then4: markedSeeds← markedSeeds + seed5: if there is at least one < 7 mm distant seed then6: markedSeeds← markedSeeds + seed7:8: procedure seed correction9: r← default neighborhood size10: for a randomly chosen marked seed in markedSeeds do11: search:12: availableSpots← find all spots on Ψ in a neighborhood of r13: for all spots do14: if spot is a forbidden spot on the grid then15: availableSpots← availableSpots - spot16: if min distance to the closest seed < 7 mm then17: availableSpots← availableSpots - spot18: if no needle exist in the new spot then19: availableSpots← availableSpots - spot20:21: if availableSpots is empty then22: r ← r + 1.23: goto search.24: else25: seed← update with a randomly chosen spot in availableSpots865.3. Experiments and Results5.3 Experiments and ResultsA cohort of 590 patients who underwent prostate brachytherapy at the Van-couver Cancer Centre (VCC), is used to evaluate the proposed automaticpreplanning framework. For each patient, PTV contours and actual seedarrangement used for the procedure are collected. All patients were treatedwith a plan based on the standard monotherapy dose prescription of 144 GymPD.We performed our experiments by dividing the cohort of patients into80% training and 20% test cases. A 5-fold cross validation approach wasimplemented to measure performance of the proposed framework. Objectivefunction coefficients were heuristically chosen as α = 0.6, β = 0.7, γ = 1.0and δ = 0.1. The rationale behind choosing importance weights is to treatV200 as the most significant term and reduce the coefficients for V150 and V100relatively. This is due to the fact that V200 is closer to the geometrical seedsarrangement in contrast to V100 which represents overall dose coverage. Weassign lower weight to the term representing the total number of needles, topush the optimizer towards seed reconfiguration rather than minimizationof total needle count. Overall performance of the optimizer is dependent onproper weight assignment. Heuristically chosen values produce sub-optimalresults.The performance analysis consists of two parts: 1) plan estimation fol-lowed by optimization, and 2) optimization of randomly generated plans. Weuse the same parameters for the optimizer in the two experiments. Totalnumber of required seeds, σ, is determined from the PTV size as indicatedin Sec. 5.2.2.We repeat our experiments for 95% confidence interval extremes, σ− andσ+, to determine the number of required seeds in both analyses; hence, wereport performance of the plans for all three runs (i.e., with σ, σ− and σ+) foreach test case. To tune the hyper-parameters of the proposed joint sparsedictionary model, we randomly partitioned the cohort into 80% trainingand 20% validation datasets. An exhaustive grid search on dictionary sizeand level of sparsity resulted in optimal number of 20 atoms with 5 sparsecoefficients.Figure 5.6 shows general statistics of the plan quality indicators V100,V150 and V200 obtained from the optimization process using estimated planand randomly generated plan. Hatched area shows the recommended zonefor each parameter according to clinical guidelines. Results tabulated inTable 5.2 further show performance of the proposed optimization algorithmin terms of dosimetry parameters and total number of needles. It is seen that875.4. Discussion and Conclusionthe optimizer can further fine-tune the actual plans towards an improved setof quality indicators.Table 5.2: Mean and standard deviation of the plan quality metrics calcu-lated for actual, estimated and random plans with and without optimization.V100(%) V150(%) V200(%) D90(Gy) #NeedlesNo Optimization:Actual 96.39±1.22 52.60±3.78 19.35±4.40 167.99±3.14 25±3Proposed Estimator 91.71±5.56 54.99±7.57 25.48±5.85 156.27±14.89 26±4Randomly Initialized 48.50±15.78 16.00±6.76 7.78±3.03 91.45±18.75 59±5With Optimization:Actual 97.13±1.35 52.15±2.91 17.80±2.25 168.39±3.82 26±3Proposed Estimator 97.12±1.95 52.32±3.48 18.03±2.54 167.79±6.62 28±4Randomly Initialized 95.04±2.60 48.94±4.45 17.46±2.76 159.93±7.12 31±5A successful plan needs to meet all or most clinical standards; therefore,we compare the percentage of successful preplans in terms of recommendedclinical guidelines and actual delivered plans. Our observation shows 47%and 14% success rate for the optimizer initialized by estimated plan andrandom plan, respectively. These rates rise up to 86% and 65% when bestof the optimization result is chosen for σ, σ− and σ+.We conducted a paired t-test analysis with the null hypothesis indicatingwhether calculated dosimetry metrics of the plans obtained from differentapproaches, i.e. learning-based and randomly initialized seed plans, havethe same mean values. We observed a statically significant difference be-tween obtained quality indices with p-value< 0.001 indicating that the nullhypothesis is rejected. Since mean values for all quality indices for ran-domly initialized seed plans were below the clinical guidelines, we obtaineda significant improvement by initializing seed plans from the learning-basedalgorithm.According to the results shown in Table 5.2 the proposed estimator pro-duces plans with quality indicators within the same range of the actualdelivered plans; as paired t-test analysis for mean values of all four dosime-try parameters between actual and proposed plans shows equivalence of themean values with p-value< 0.001.5.4 Discussion and ConclusionIn this work, an automatic prostate brachytherapy treatment planning frame-work is introduced that incorporates previous successful treatment knowl-edge in a model by capturing joint sparse patterns between PTV and seed885.4. Discussion and Conclusion 020406080100% V100 V150 V200Actual Plan-+(a) 020406080100% V100 V150 V200Actual Plan-+(b)Figure 5.6: Comparison between statistics of plan quality metrics for a)estimated plan followed by optimization, and b) randomly generated planfollowed by optimization. Number of seeds used in both methods are ob-tained from polynomial regression model.895.4. Discussion and Conclusionconfiguration. Results of the model, i.e. plan estimator, is further optimizedby a stochastic optimizer that uses simulated annealing architecture and anovel cost function to apply clinical guidelines and improve the plan qualityindicators. Initialization of the optimizer from the estimated plan reducesthe very large search space for all possible plans to a neighborhood of theinitially estimated plan.The proposed system is evaluated on a retrospective dataset of treat-ment records at the Vancouver Cancer Centre. We evaluate performance ofthe system by reporting the statistics of the plan quality indicators againstthe actual delivered plan, and show a very high success rate to converge toacceptable plans in terms of clinical guidelines. We show that the automati-cally generated plans require same range of needles to deliver the plan whencompared with the actual plans used for treatment.Each of the two main components of the proposed framework are at-tributed to some limitations. In this study, learning-based sparse dictionarygeneration is applied to learn geometrical distribution of the seeds irrespec-tive to their activation and dose distribution since our training dataset al-lowed. It is more desired to introduce a mathematical model which is trainedfor dose distribution patterns, otherwise multiple models for different seedactivations may be needed. The proposed optimization algorithm tunes theseed arrangements according to the strict clinical guidelines; however, in-troduces a large set of parameters to the proposed preplanning framework.Optimizer parameter tuning would benefit from sensitivity analysis.Our proposed method can potentially aid the clinicians in training, de-cision support or suggesting plans which can be further modified. The pro-posed framework can be used as a tool to improve consistency across variouscentres and clinicians by automatically providing a reference plan that is incompliance with knowledge existed in one unified training dataset. Theoptimization algorithm can also be used to reduce the number of needlesrequired to deliver the plan.The framework is computationally efficient, currently running on a com-modity personal computer with an unoptimized MATLABR©code within2 minutes. Hence, it can be further optimized for intra-operative applica-tions. The framework can also adapt to different institutions, requirementsor datasets, and is capable of predicting a clinically acceptable plan thatcaptures an implicit knowledge of a clinician or a group of clinicians, so thatclinician-specific planning can be performed.90Chapter 6Conclusion and Future WorkCurrent flow of the information for prostate brachytherapy preplanning pro-cedure is prone to subjective errors, due to poor visibility of the targetanatomy in TRUS images and high dependency to expertise of the physi-cian. Accumulation of the errors in preplanning processes can potentiallyinfluence treatment outcome and rate of morbidity. In this thesis, in aneffort to reduce the need for user interaction, we have introduced methodsand algorithms to automate sequential processes involved in brachytherapypreplanning. Various information fusion algorithms have been introduced torepresent joint relation between volumetric information elements of preplan-ning in the form of mathematical and statistical models. These models aretrained based on previous treatment records obtained from the VancouverCancer Centre that have shown to result in an excellent treatment outcome.In Chapter 2, we introduced a CTV delineation algorithm based on fu-sion of multiple atlases. Prior knowledge in delineation of CTV contoursis incorporated in a set of atlases and is propagated through a non-rigiddeformable registration algorithm to a new patient’s TRUS volume. Theproposed method generated clinically acceptable results when compared tothe gold standards; however, it lacks the requirements for real-time appli-cations such as intra-operative dosimetry and plan correction, because ofthe computation expenses. In search for a fast and reliable fusion algo-rithm, we formulated the segmentation problem with a linear joint modelthat encapsulates the joint patterns between TRUS intensity volumes andthe corresponding contours in a set of joint components. In Chapter 3, weadapted the independent component analysis to estimate joint componentsand presented its application in estimation of CTV, PTV and mPD con-tours. In Chapter 4, we extended the proposed methodology by introducingsparsity to the modeling process. Further, we proposed a framework forsimultaneous estimations of two planning labels, i.e. CTV and PTV. Theproposed approach has been evaluated on a large dataset of 590 treatmentcase.In Chapter 5, we proposed a fusion algorithm for estimation of the seedarrangement. We introduced a two-part framework that firstly generates91Chapter 6. Conclusion and Future Workinitial seed distribution utilizing the joint sparse model, and secondly re-arranges seeds using a greedy optimization algorithm. Generated plans arecompared against actual plans in terms of plan quality indicators. We con-ducted adaptability analysis to investigate a feature of the proposed method-ology in capturing dataset-specific knowledge. We concluded that the pro-posed framework can generate better plans in terms of dosimetry parameterswhen trained and evaluated on a dataset of the same physician.In conclusion, in this thesis, we aimed to automate various componentsof the preplanning procedure for prostate brachytherapy. We investigatedvarious atlas-based and join statistical modeling approaches, while incor-porating clinical guidelines. We achieved our goal by analyzing a largeretrospective brachytherapy treatment planning dataset at the BC CancerAgency. We demonstrated that real-time estimation of CTV and PTV isfeasible. Further, we solved the inverse problem of seed placement by us-ing a learning-based initial seed distribution to obtain a plan that closelycomplies with the clinical guidelines.The contributions of this thesis are summarized as follows:• A multi-atlas fusion framework is introduced for automatic delineationof the CTV from TRUS images. A dataset of a priori segmented ul-trasound images, i.e. atlases, is registered to a target image. Cor-responding labels of the atlases are fused based on pairwise shapesimilarity. Evaluation of the proposed segmentation approach on a setof transrectal prostate volume studies produces segmentation resultsthat are within the range of observer variability when compared to asemi-automatic segmentation technique that is routinely used at theVCC.• A fusion criteria is introduced to sort and rank the atlases in themulti-atlas-based segmentation approach. We propose pairwise atlasagreement factor that combines an image-similarity metric and simi-larity between a priori segmented contours. This factor is used in anatlas selection algorithm to prune the dataset before combining theatlas contours to produce a consensus segmentation.• A fast and efficient fusion framework based on ICA is introduced forestimation of single or multiple planning contours, i.e. the CTV, thePTV and the mPD. Generally, this linear decomposition algorithm istrained on a set of complete joint observation vectors, and is tested on asingle partial observation vector to estimate missing parts of the obser-vation vector. We introduce a method to calculate a global threshold926.1. Future Workvalue from training dataset that converts the estimated probabilitymaps into final binary volumes.• A sparsity-based fusion framework is introduced for simultaneous esti-mation of the CTV and the PTV. This framework is a generalizationof the ICA-based modeling approach that is presented to automate de-lineation of CTV and PTV by finding sparse joint patterns extractedfrom previous treatment records.• An optimization algorithm is proposed to produce seed arrangementw.r.t. the clinical guidelines at the VCC. This optimizer is designedon top of a simulated annealing architecture and encodes some keyguidelines in both objective function and the annealing process.• We propose a learning-based seed planning framework that combinesthe proposed optimizer with the sparsity-based fusion framework. Thelearning-based fusion framework provides seed configuration based onsparse joint patterns in a planning dataset. The initial seed arrange-ment is further perturbed according to the clinical guidelines using thein-house optimizer. This plan estimation framework has been evalu-ated on a large dataset of treatment records collected from the VCCand shows acceptable results based on clinically suggested plan qualityindicators.6.1 Future WorkNovel methods have been presented in this thesis for automation of theprostate brachytherapy planning procedure using information fusion andmachine learning techniques. A number of interesting areas of research canbe suggested as follows:• The proposed CTV segmentation using multiple atlases is essentiallybuilt on top of an intensity-based registration process. Image reg-istration is more challenging in the context of TRUS images mainlybecause of the poor visibility of the CTV near the base and the apex.Definition of the CTV is also a subjective concept which limits theproposed approach to correctly estimate the target contours for a newcase. To alleviate this, a very large dataset is required that includesvariety of CTV variations. hence we recommend:936.1. Future Work– Using a larger dataset of atlases. This increases likelihood ofhaving more similar cases in the atlas dataset for a new patientTRUS volume.– To regularize intensity-based deformable registration algorithmbased on the location proximity to the base or the apex.– In case of a very large dataset of atlases, the proposed frameworkwould benefit from pre-registration grouping of atlases w.r.t. asimilarity metric. This can be implemented either generally us-ing pair- or group-wise clustering algorithms, or particularly byincluding new patient’s TRUS volume in similarity calculation.– Our CPU-based non-optimized implementation of the frameworkproduces segmentation results in few minutes. Since atlas regis-tration processes are performed individually and independently,our implementation can potentially be improved by orders ofmagnitude if it utilizes graphical processing unit (GPU) for theregistration algorithm. An improved implementation can poten-tially be used in real-time segmentation applications such as in-traoperative dosimetry and prostate guided biopsy.• Our proposed joint CTV and PTV estimation framework takes ad-vantage of sparsity constraints in generating joint fusion-based mod-els. The method is efficient, fast and reliable; however, accuracy ofestimation can be improved by considering the following points:– A larger dataset and a more generic model that can be appliedto different size and shape prostate cases. Bootstrapping anddataset pruning based on some clinical labels such as size, shape,biomarkers or image quality provides one or multiple joint modelsthat represent clinically appealing variations of the patients.– We propose a method to convert estimated probability maps ofthe target volumes into final binary estimation of the target vol-umes. This approach uses a dataset (training) dependent globalthreshold value for all the query TRUS volumes. We recommendusing a case-specific threshold value that can use priors from thetraining dataset to estimate optimum threshold value per eachquery case.– The sparsity-based fusion algorithm assumes all the volumetricinformation representation of the context are roughly pre-aligned.This is guaranteed in case of the VCC dataset, since clinical guide-lines oblige the physician to maintain the prostate in the middle946.1. Future Workof the TRUS images while moving the probe axially from the baseto the apex. We suggest to add a pre-processing step that per-forms a rigid alignment of the training and the test cases to ensurethis assumption is held true independent of the data acquisitionprotocol.• The proposed target volume segmentation approaches have to be inte-grated into the current clinical workflow, e.g. at the VCC. It is recom-mended to provide an efficient implementation of the algorithms thatcan communicate with the current software in clinics with minimumchange in the existing software infrastructure.• A definition of confidence for segmentation is highly desirable. Thisis mostly beneficial for clinical integration, as it provides a confidencelevel for contours to help clinicians to perform corrections. Confidencecan be defined from level of similarity to those cases included in thetraining set.• Methods proposed and reviewed in this thesis are based on a linearmodeling approach, while complexity of the context is captured viaextracted joint patterns. Utilizing some non-linear methods may im-prove accuracy of the model in capturing the joint patterns in thefused information space.• In the course of investigating fusion algorithms, we see ICA and sparseanalysis as special cases of a general modeling approach. We recom-mend to review application of some unsupervised techniques that en-code observation vectors into natural clusters such as restricted Bultz-man machines (RBM) [98], auto-associative networks and random for-est [99].• Although adaptability can be one of the key features of the proposedmethodology, a remaining question to be addressed as future workis to determine how applicable this methodology is to other datasetsoutside VCC.– The dataset used in this study is limited to a large cohort ofpatients imaged by the same brand ultrasound machine whichhas been widely used across BC. However, any other dataset thatprovides an extent of the anatomy which observations are reason-ably standardized can be used for joint pattern analysis. Imaging956.1. Future Workparameters and hardware related artifacts might appear as un-wanted complex patterns that may influence performance of thesegmentation. It is clear that removing those disruptive patternsby standardizing the imaging protocol or splitting the trainingdataset provides better resources for training models with morepowerful joint patterns expressing the segmentation problem.– Joint patterns learned by the model represent the segmentationrelated information in the training dataset. Therefore, a newpatient with incompatible imaging data (e.g. non-axial view, dif-ferent spacing or unusual imaging parameters) can be detectedfrom reconstruction error after projection onto the joint compo-nents. Eventually, an inclusion criteria can be defined to prunethe outliers in test time.• Plan estimation approach presented in this thesis shows feasibility ofautomatic seed planning using a two component automation frame-work: 1) A fusion-based learning component that estimated a seedconfiguration in accordance with the priors in the training dataset;and 2) an optimization algorithm that strictly enforces the clinicalguidelines per case. Current approach limits the optimizer to searchfor a clinically acceptable seed arrangement in the neighborhood of theestimated seed configuration. This approach can be pushed more to-wards the priors by mixing the learning-based model and optimizationalgorithm. We recommend to guide the optimization process, per it-eration, to follow patterns more compliant with the joint model. Thiscan be implemented in the annealing process, where the seed arrange-ment is modified slightly towards the clinical guidelines. We suggestto introduce more optimization constraints to the objective functionaccording to the VCC guidelines:– Total number of seeds outside PTV is desired to be minimum.Hence, priors of such factor can be used in the objective functionto minimize this value. Our presented framework generates planswith greater number of seeds outside PTV for the larger sizeprostates, while this value is less for the smaller size prostateswhen compared to the actual plans. A new weighted term can beintroduced to the objective function to apply this constraint.– The proposed approach generates plans with higher entropy interms of seed arrangement when compared to the actual plans.This is observable when analyzing plans in the coronal views,966.1. Future Workwhich show less regularity in the peripheral zone. For example,needles with single seeds are not common in actual plans; how-ever, the optimizer is not constrained to avoid use of such needles.We suggest to use a grid template probability map that provides aguide for placement of regularly loaded strands or special loadedneedles.– It is not usually recommended to place three vertical needles ina column in anterior regions. Although limit of three verticalneedles has been introduced to the optimizer, it is not constrainedlocation-wise.– One of the reasons actual plans express more regularity comparedwith the proposed method is the fact that actual preplanning pro-cess starts off by placing regular strands in desired grid locations.Hence, we suggest to generate a separate model for needle con-figuration based on PTV and utilize it in the annealing processor in initializing seed arrangement.97Bibliography[1] S. S. Mahdavi, N. Chng, I. Spadinger, W. J. Morris, and S. E. Salcud-ean, “Semi-automatic segmentation for prostate interventions,” Med.Image Anal., vol. 15, no. 2, pp. 226–237, 2011.[2] C. K. e. Howlader N, Noone AM, Krapcho M, Garshell J, Neyman N,Altekruse SF, Kosary CL, Yu M, Ruhl J, Tatalovich Z, Cho H, MariottoA, Lewis DR, Chen HS, Feuer EJ, “SEER Cancer Statistics Review,1975-2010,” Bethesda, MD, 2012.[3] Advisory Committee on Cancer Statistics, “Canadian Cancer Statis-tics,” 2013.[4] K. A. Hinnen, J. J. Battermann, J. G. H. van Roermund, M. A. Moer-land, I. M. Ju¨rgenliemk-Schulz, S. J. Frank, and M. van Vulpen, “Long-term biochemical and survival outcome of 921 patients treated withI-125 permanent prostate brachytherapy,” Int. J. Radiat. Oncol. Biol.Phys., vol. 76, no. 5, pp. 1433–8, Apr. 2010.[5] W. J. Morris, M. Keyes, I. Spadinger, W. Kwan, M. Liu, M. McKen-zie, H. Pai, T. Pickles, and S. Tyldesley, “Population-based 10-yearoncologic outcomes after low-dose-rate brachytherapy for low-risk andintermediate-risk prostate cancer,” Cancer, vol. 119, no. 8, pp. 1537–46,Apr. 2013.[6] M. Keyes, W. J. Morris, I. Spadinger, C. Araujo, A. Cheung, N. Chng,J. Crook, R. Halperin, V. Lapointe, S. Miller, H. Pai, and T. Pick-les, “Radiation oncology and medical physicists quality assurance inBritish Columbia Cancer Agency Provincial Prostate BrachytherapyProgram,” Brachytherapy, vol. 12, no. 4, pp. 343–55, 2013.[7] J. E. Sylvester, P. D. Grimm, S. M. Eulau, R. K. Takamiya, andD. Naidoo, “Permanent prostate brachytherapy preplanned technique:The modern Seattle method step-by-step and dosimetric outcomes,”Brachytherapy, vol. 8, no. 2, pp. 197–206, 2009.98Bibliography[8] Y. Yu, L. L. Anderson, Z. Li, D. E. Mellenberg, R. Nath, M. C. Schell,F. M. Waterman, A. Wu, and J. C. Blasko, “Permanent prostate seedimplant brachytherapy: report of the American Association of Physi-cists in Medicine Task Group No. 64,” Med. Phys., vol. 26, no. 10, pp.2054–76, Oct. 1999.[9] W. L. Smith, C. Lewis, G. Bauman, G. Rodrigues, D. D’Souza, R. Ash,D. Ho, V. Venkatesan, D. Downey, and A. Fenster, “Prostate volumecontouring: a 3D analysis of segmentation using 3DTRUS, CT, andMR,” Int. J. Radiat. Oncol. Biol. Phys., vol. 67, no. 4, pp. 1238–47,Mar. 2007.[10] S. Nag, W. Bice, K. DeWyngaert, B. Prestidge, R. Stock, and Y. Yu,“The American Brachytherapy Society recommendations for perma-nent prostate brachytherapy postimplant dosimetric analysis,” Int. J.Radiat. Oncol. Biol. Phys., vol. 46, no. 1, pp. 221–30, Jan. 2000.[11] Z. Chen, N. Yue, X. Wang, K. B. Roberts, R. Peschel, and R. Nath,“Dosimetric effects of edema in permanent prostate seed implants: arigorous solution,” Int. J. Radiat. Oncol. Biol. Phys., vol. 47, no. 5, pp.1405–19, Jul. 2000.[12] N. Yue, Z. Chen, R. Peschel, A. P. Dicker, F. M. Waterman, andR. Nath, “Optimum timing for image-based dose evaluation of 125Iand 103PD prostate seed implants,” Int. J. Radiat. Oncol. Biol. Phys.,vol. 45, no. 4, pp. 1063–72, Nov. 1999.[13] M. Keyes, C. Macaulay, M. Hayes, J. Korbelik, W. J. Morris, andB. Palcic, “DNA ploidy measured on archived pretreatment biopsymaterial may correlate with prostate-specific antigen recurrence afterprostate brachytherapy,” Int. J. Radiat. Oncol. Biol. Phys., vol. 86,no. 5, pp. 829–34, Aug. 2013.[14] W. D. Richard and C. G. Keen, “Automated texture-based segmen-tation of ultrasound images of the prostate,” Comput. Med. ImagingGraph., vol. 20, no. 3, pp. 131–40, 1996.[15] S. D. Pathak, V. Chalana, D. R. Haynor, and Y. Kim, “Edge-guidedboundary delineation in prostate ultrasound images,” IEEE Trans.Med. Imaging, vol. 19, no. 12, pp. 1211–9, Dec. 2000.99Bibliography[16] H. M. Ladak, F. Mao, Y. Wang, D. B. Downey, D. A. Steinman, andA. Fenster, “Prostate boundary segmentation from 2D ultrasound im-ages,” Med. Phys., vol. 27, pp. 1777–88, 2000.[17] P. Abolmaesumi and M. R. Sirouspour, “An interacting multiple modelprobabilistic data association filter for cavity boundary extraction fromultrasound images,” IEEE Trans. Med. Imaging, vol. 23, no. 6, pp.772–784, 2004.[18] Y. Wang, H. N. Cardinal, D. B. Downey, and A. Fenster, “Semi-automatic three-dimensional segmentation of the prostate using two-dimensional ultrasound images,” Med. Phys., vol. 30, pp. 887–97, 2003.[19] M. Ding, B. Chiu, I. Gyacskov, X. Yuan, M. Drangova, D. B. Downey,and A. Fenster, “Fast prostate segmentation in 3D TRUS images basedon continuity constraint using an autoregressive model,” Med. Phys.,vol. 34, pp. 4109–25, 2007.[20] C. Knoll, M. Alcan˜iz, V. Grau, C. Monserrat, and M. C. Juan, “Out-lining of the prostate using snakes with shape restrictions based on thewavelet transform,” Pattern Recognit., vol. 32, no. 10, pp. 1767–1781,1999.[21] W. Qiu, J. Yuan, E. Ukwatta, D. Tessier, and A. Fenster, “Rotational-slice-based prostate segmentation using level set with shape constraintfor 3D end-firing TRUS guided biopsy,” Med. Image Comput. Comput.Interv. MICCAI, vol. 15, no. Pt 1, pp. 537–44, Jan. 2012.[22] A. C. Hodge, A. Fenster, D. B. Downey, and H. M. Ladak, “Prostateboundary segmentation from ultrasound images using 2D active shapemodels: Optimisation and extension to 3D,” Comput. Methods Pro-grams Biomed., vol. 84, no. 2, pp. 99–113, 2006.[23] L. Gong, S. D. Pathak, D. R. Haynor, P. S. Cho, and Y. Kim, “Para-metric shape modeling using deformable superellipses for prostate seg-mentation,” IEEE Trans. Med. Imaging, vol. 23, no. 3, pp. 340–349,2004.[24] S. S. Mahdavi, I. Spadinger, N. Chng, S. E. Salcudean, and W. J. Mor-ris, “Semiautomatic segmentation for prostate brachytherapy: dosimet-ric evaluation,” Brachytherapy, vol. 12, no. 1, pp. 65–76, Jan. 2013.100Bibliography[25] S. S. Mahdavi, M. Moradi, W. J. Morris, S. L. Goldenberg, and S. E.Salcudean, “Fusion of ultrasound B-mode and vibro-elastography im-ages for automatic 3D segmentation of the prostate,” IEEE Trans. Med.Imaging, vol. 31, no. 11, pp. 2073–2082, Jul. 2012.[26] S. Ghose, A. Oliver, J. Mitra, R. Mart´ı, X. Llado´, J. Freixenet,D. Sidibe´, J. C. Vilanova, J. Comet, and F. Meriaudeau, “A super-vised learning framework of statistical shape and probability priors forautomatic prostate segmentation in ultrasound images,” Med. ImageAnal., vol. 17, no. 6, pp. 587–600, Aug. 2013.[27] Y. Zhan and D. Shen, “Deformable segmentation of 3D ultrasoundprostate images using statistical texture matching method,” IEEETrans. Med. Imaging, vol. 25, no. 3, pp. 256–272, 2006.[28] P. Yan, S. Xu, B. Turkbey, and J. Kruecker, “Discrete DeformableModel Guided by Partial Active Shape Model for TRUS Image Seg-mentation,” IEEE Trans. Biomed. Eng., vol. 57, no. 5, pp. 1158–66,May 2010.[29] W. Qiu, J. Yuan, E. Ukwatta, and A. Fenster, “Rotationally resliced3D prostate TRUS segmentation using convex optimization with shapepriors,” Med. Phys., vol. 42, no. 2, pp. 877–91, Feb. 2015.[30] P. Wu, Y. Liu, Y. Li, and B. Liu, “Robust Prostate Segmentation Us-ing Intrinsic Properties of TRUS Images,” IEEE Trans. Med. Imaging,vol. 34, no. 6, pp. 1321–1335, Jan. 2015.[31] Q. Xie and D. Ruan, “Low-complexity atlas-based prostate segmen-tation by combining global, regional, and local metrics,” Med. Phys.,vol. 41, no. 4, p. 041909, Apr. 2014.[32] H. Akbari, X. Yang, L. V. Halig, and B. Fei, “3D segmentation ofprostate ultrasound images using wavelet transform,” in Proc. SPIE,vol. 7962, 2011, p. 79622K.[33] X. Yang and B. Fei, “3D prostate segmentation of ultrasound imagescombining longitudinal image registration and machine learning,” inProc. SPIE, vol. 8316, 2012, p. 83162O.[34] T. Heimann, M. Baumhauer, T. Simpfendo¨rfer, H. P. Meinzer, andI. Wolf, “Prostate segmentation from 3D transrectal ultrasound us-ing statistical shape models and various appearance models,” in Proc.SPIE, vol. 6914, 2008, p. 69141P.101Bibliography[35] R. Toth, J. Ribault, J. Gentile, D. Sperling, and A. Madabhushi, “Si-multaneous Segmentation of Prostatic Zones Using Active AppearanceModels With Multiple Coupled Levelsets,” Comput. Vis. Image Un-derst., vol. 117, no. 9, pp. 1051–1060, Sep. 2013.[36] J. Yang and J. S. Duncan, “3D image segmentation of deformable ob-jects with joint shape-intensity prior models using level sets,” Med.Image Anal., vol. 8, no. 3, pp. 285–94, Sep. 2004.[37] S. Klein, U. A. van der Heide, I. M. Lips, M. van Vulpen, M. Staring,and J. P. W. Pluim, “Automatic segmentation of the prostate in 3D MRimages by atlas matching using localized mutual information,” Med.Phys., vol. 35, p. 1407, 2008.[38] J. A. Dowling, J. Fripp, S. Chandra, J. P. W. Pluim, J. Lambert,J. Parker, J. Denham, P. B. Greer, and O. Salvado, “Fast automaticmulti-atlas segmentation of the prostate from 3D MR images,” ProstateCancer Imaging Image Anal. Image-Guided Interv., pp. 10–21, 2011.[39] T. Rohlfing, R. Brandt, R. Menzel, and C. R. Maurer, “Evaluationof atlas selection strategies for atlas-based image segmentation withapplication to confocal microscopy images of bee brains,” Neuroimage,vol. 21, no. 4, pp. 1428–1442, 2004.[40] C. Lu, S. Chelikani, X. Papademetris, J. P. Knisely, M. F. Milosevic,Z. Chen, D. A. Jaffray, L. H. Staib, and J. S. Duncan, “An integratedapproach to segmentation and nonrigid registration for application inimage-guided pelvic radiotherapy,” Med. Image Anal., vol. 15, no. 5, p.772, 2011.[41] O. Acosta, A. Simon, F. Monge, F. Commandeur, C. Bassirou, G. Ca-zoulat, R. de Crevoisier, and P. Haigron, “Evaluation of multi-atlas-based segmentation of CT scans in prostate cancer radiotherapy,” inBiomed. Imaging (ISBI), IEEE Int. Symp. IEEE, 2011, pp. 1966–1969.[42] X. Artaechevarria, A. Munoz-Barrutia, and C. Ortiz-de Solorzano,“Combination strategies in multi-atlas image segmentation: applica-tion to brain MR data,” IEEE Trans. Med. Imaging, vol. 28, no. 8, pp.1266–77, Aug. 2009.[43] R. A. Heckemann, J. V. Hajnal, P. Aljabar, D. Rueckert, and A. Ham-mers, “Automatic anatomical brain MRI segmentation combining label102Bibliographypropagation and decision fusion,” Neuroimage, vol. 33, no. 1, pp. 115–26, Oct. 2006.[44] P. Aljabar, R. A. Heckemann, A. Hammers, J. V. Hajnal, and D. Rueck-ert, “Multi-atlas based segmentation of brain images: atlas selectionand its effect on accuracy,” Neuroimage, vol. 46, no. 3, pp. 726–38, Jul.2009.[45] A. Redpath, “Automatic determination of needle and source positionsfor brachytherapy of the prostate using 125Iodine Rapid Strand,” Ra-diother. Oncol., vol. 64, no. 2, pp. 215–27, Aug. 2002.[46] G. Ferrari, Y. Kazareski, F. Laca, and C. E. Testuri, “A model forprostate brachytherapy planning with sources and needles position op-timization,” Oper. Res. Heal. Care, vol. 3, no. 1, pp. 31–39, Mar. 2014.[47] J. Pouliot, D. Tremblay, J. Roy, and S. Filice, “Optimization of per-manent 125I prostate implants using fast simulated annealing,” Int. J.Radiat. Oncol. Biol. Phys., vol. 36, no. 3, pp. 711–20, Oct. 1996.[48] R. R. Meyer, W. D. D’Souza, M. C. Ferris, and B. R. Thomadsen, “MIPModels and BB Strategies in Brachytherapy Treatment Optimization,”J. Glob. Optim., vol. 25, no. 1, pp. 23–42, Jan. 2003.[49] E. Lessard, S. L. S. Kwa, B. Pickett, M. Roach, and J. Pouliot, “Classsolution for inversely planned permanent prostate implants to mimican experienced dosimetrist,” Med. Phys., vol. 33, no. 8, p. 2773, Jul.2006.[50] Y. Yu, J. B. Zhang, R. A. Brasacchio, P. G. Okunieff, D. J. Rubens,J. G. Strang, A. Soni, and E. M. Messing, “Automated treatment plan-ning engine for prostate seed implant brachytherapy,” Int. J. Radiat.Oncol. Biol. Phys., vol. 43, no. 3, pp. 647–52, Feb. 1999.[51] E. K. Lee, R. J. Gallagher, D. Silvern, C. S. Wuu, and M. Zaider,“Treatment planning for brachytherapy: an integer programmingmodel, two computational approaches and experiments with perma-nent prostate implant planning,” Phys. Med. Biol., vol. 44, no. 1, pp.145–65, Jan. 1999.[52] G. Yang, L. E. Reinstein, S. Pai, Z. Xu, and D. L. Carroll, “A newgenetic algorithm technique in optimization of permanent 125I prostateimplants,” Med. Phys., vol. 25, no. 12, pp. 2308–15, Dec. 1998.103Bibliography[53] R. J. Gallagher and E. K. Lee, “Mixed integer programming optimiza-tion models for brachytherapy treatment planning,” Proc. AMIA Annu.Fall Symp., pp. 278–82, Jan. 1997.[54] S. Nouranian, S. S. Mahdavi, I. Spadinger, W. Morris, S. Salcudean,and P. Abolmaesumi, “A Multi-Atlas-Based Segmentation Frameworkfor Prostate Brachytherapy,” IEEE Trans. Med. Imaging, vol. 34, no. 4,pp. 950–961, Dec. 2015.[55] S. Liao, Y. Gao, J. Lian, and D. Shen, “Sparse patch-based label prop-agation for accurate prostate localization in CT images,” IEEE Trans.Med. Imaging, vol. 32, no. 2, pp. 419–34, Feb. 2013.[56] M. R. Sabuncu, B. T. T. Yeo, K. Van Leemput, B. Fischl, and P. Gol-land, “A generative model for image segmentation based on label fu-sion,” IEEE Trans. Med. Imaging, vol. 29, no. 10, pp. 1714–29, Oct.2010.[57] O. Commowick, A. Akhondi-Asl, and S. K. Warfield, “Estimating a ref-erence standard segmentation with spatially varying performance pa-rameters: local MAP STAPLE,” IEEE Trans. Med. Imaging, vol. 31,no. 8, pp. 1593–606, Aug. 2012.[58] A. Akhondi-Asl, L. Hoyte, M. Lockhart, and S. Warfield, “A Loga-rithmic Opinion Pool Based STAPLE Algorithm for the Fusion of Seg-mentations With Associated Reliability Weights,” IEEE Trans. Med.Imaging, Jun. 2014.[59] T. R. Langerak, F. F. Berendsen, U. A. Van der Heide, A. N. T. J.Kotte, and J. P. W. Pluim, “Multiatlas-based segmentation with pre-registration atlas selection,” Med. Phys., vol. 40, no. 9, p. 091701, Sep.2013.[60] S. Nouranian, S. S. Mahdavi, I. Spadinger, W. J. Morris, S. E. Salcud-ean, and P. Abolmaesumi, “An Automatic Multi-atlas Segmentationof the Prostate in Transrectal Ultrasound Images Using Pairwise AtlasShape Similarity,” in Med. Image Comput. Comput. Interv. MICCAI,2013, vol. 8150, pp. 173–180.[61] ——, “Multi-atlas-based automatic 3D segmentation for prostatebrachytherapy in transrectal ultrasound images,” Proc. SPIE Med.Imaging, vol. 8671, pp. 86 710O–86 710O–7, 2013.104Bibliography[62] J. P. Thirion, “Image matching as a diffusion process: an analogy withMaxwell’s demons,” Med. Image Anal., vol. 2, no. 3, pp. 243–260, 1998.[63] X. Pennec, P. Cachier, and N. Ayache, “Understanding the demons al-gorithm: 3D non-rigid registration by gradient descent,” in Med. ImageComput. Comput. Interv. MICCAI. Springer, 1999, pp. 597–605.[64] T. Vercauteren, X. Pennec, A. Perchant, and N. Ayache, “Non-parametric Diffeomorphic Image Registration with Demons Algo-rithm,” in Med. Image Comput. Comput. Interv., ser. Lecture Notesin Computer Science, N. Ayache, S. Ourselin, and A. J. Maeder, Eds.,vol. 10, no. Pt 2, Asclepios Research Group, INRIA Sophia-Antipolis,France. Springer, 2007, pp. 319–326.[65] S. Khallaghi, C. G. M. Leung, K. Hastrudi-Zaad, P. Foroughi,C. Nguan, and P. Abolmaesumi, “Experimental validation of an intra-subject elastic registration algorithm for dynamic-3D ultrasound im-ages,” Med. Phys., vol. 39, no. 9, pp. 5488–97, Sep. 2012.[66] T. De Silva, A. Fenster, J. Bax, C. Romagnoli, J. Izawa, J. Samara-bandu, and A. D. Ward, “Quantification of prostate deformation due toneedle insertion during TRUS-guided biopsy: comparison of hand-heldand mechanically stabilized systems,” Med. Phys., vol. 38, no. 3, pp.1718–31, Mar. 2011.[67] I. Isgum, M. Staring, A. Rutten, M. Prokop, M. A. Viergever, andB. van Ginneken, “Multi-atlas-based segmentation with local decisionfusion–application to cardiac and aortic segmentation in CT scans,”IEEE Trans. Med. Imaging, vol. 28, no. 7, pp. 1000–10, Jul. 2009.[68] T. R. Langerak, U. A. van der Heide, A. N. T. J. Kotte, M. A. Viergever,M. van Vulpen, and J. P. W. Pluim, “Label fusion in atlas-based seg-mentation using a selective and iterative method for performance levelestimation (SIMPLE),” IEEE Trans. Med. Imaging, vol. 29, no. 12, pp.2000–8, Dec. 2010.[69] S. K. Warfield, K. H. Zou, and W. M. Wells, “Simultaneous truth andperformance level estimation (STAPLE): an algorithm for the valida-tion of image segmentation,” IEEE Trans. Med. Imaging, vol. 23, no. 7,pp. 903–921, 2004.[70] T. Rohlfing, D. B. Russakoff, and C. R. Maurer, “Performance-based classifier combination in atlas-based image segmentation using105Bibliographyexpectation-maximization parameter estimation,” IEEE Trans. Med.Imaging, vol. 23, no. 8, pp. 983–94, Aug. 2004.[71] A. J. Asman and B. A. Landman, “Formulating spatially varying per-formance in the statistical fusion framework,” IEEE Trans. Med. Imag-ing, vol. 31, no. 6, pp. 1326–36, Jun. 2012.[72] J. Nocedal, “Updating Quasi-Newton Matrices with Limited Storage,”Math. Comput., vol. 35, no. 151, pp. 773–782, 1980.[73] T. Rohlfing, “Image similarity and tissue overlaps as surrogates forimage registration accuracy: widely used but unreliable,” IEEE Trans.Med. Imaging, vol. 31, no. 2, pp. 153–63, Feb. 2012.[74] A. J. Asman, F. W. Bryan, S. A. Smith, D. S. Reich, and B. A. Land-man, “Groupwise multi-atlas segmentation of the spinal cord’s internalstructure,” Med. Image Anal., vol. 18, no. 3, pp. 460–71, Apr. 2014.[75] I. B. Tutar, S. D. Pathak, L. Gong, P. S. Cho, K. Wallner, and Y. Kim,“Semiautomatic 3D prostate segmentation from TRUS images usingspherical harmonics,” IEEE Trans. Med. Imaging, vol. 25, no. 12, pp.1645–1654, 2006.[76] M. Modat, G. R. Ridgway, Z. A. Taylor, M. Lehmann, J. Barnes, D. J.Hawkes, N. C. Fox, and S. Ourselin, “Fast free-form deformation us-ing graphics processing units,” Comput. Methods Programs Biomed.,vol. 98, no. 3, pp. 278–284, 2010.[77] S. Nouranian, M. Ramezani, S. S. Mahdavi, I. Spadinger, W. J. Morris,S. E. Salcudean, and P. Abolmaesumi, “Fast Prostate Segmentationfor Brachytherapy based on Joint Fusion of Images and Labels,” Proc.SPIE Med. Imaging, vol. 9036, pp. 90 361A–90 361A–7, 2014.[78] V. D. Calhoun, T. Adali, and J. Liu, “A feature-based approach tocombine functional MRI, structural MRI and EEG brain imaging data,”in Eng. Med. Biol. Soc. 28th Annu. Int. Conf. IEEE, vol. 1, Jan. 2006,pp. 3672–5.[79] V. D. Calhoun, T. Adali, N. R. Giuliani, J. J. Pekar, K. A. Kiehl,and G. D. Pearlson, “Method for multimodal analysis of independentsource differences in schizophrenia: combining gray matter structuraland auditory oddball functional data,” Hum. Brain Mapp., vol. 27,no. 1, pp. 47–62, Jan. 2006.106Bibliography[80] K. Specht, R. Zahn, K. Willmes, S. Weis, C. Holtel, B. J. Krause,H. Herzog, and W. Huber, “Joint independent component analysis ofstructural and functional images reveals complex patterns of functionalreorganisation in stroke aphasia,” Neuroimage, vol. 47, no. 4, pp. 2057–63, Oct. 2009.[81] J. Liu, G. Pearlson, A. Windemuth, G. Ruano, N. I. Perrone-Bizzozero,and V. D. Calhoun, “Combining fMRI and SNP data to investigateconnections between brain function and genetics using parallel ICA,”Hum. Brain Mapp., vol. 30, no. 1, pp. 241–55, Jan. 2009.[82] S. Nouranian, M. Ramezani, S. S. Mahdavi, I. Spadinger, W. J. Morris,S. E. Salcudean, and P. Abolmaesumi, “Data fusion for planning targetvolume and isodose prediction in prostate brachytherapy,” pp. 94 151I–94 151I–7, 2015.[83] W. J. Morris, I. Spadinger, M. Keyes, J. Hamm, M. McKenzie, andT. Pickles, “Whole prostate D90 and V100: a dose-response analysis of2000 consecutive (125)I monotherapy patients,” Brachytherapy, vol. 13,no. 1, pp. 32–41, Jan. 2014.[84] M. J. Rivard, B. M. Coursey, L. A. DeWerd, W. F. Hanson, M. S. Huq,G. S. Ibbott, M. G. Mitch, R. Nath, and J. F. Williamson, “Updateof AAPM Task Group No. 43 Report: A revised AAPM protocol forbrachytherapy dose calculations,” Med. Phys., vol. 31, no. 3, pp. 633–74, Mar. 2004.[85] T. W. Lee, M. Girolami, and T. J. Sejnowski, “Independent componentanalysis using an extended infomax algorithm for mixed subgaussianand supergaussian sources,” Neural Comput., vol. 11, no. 2, pp. 417–41, Feb. 1999.[86] J. Himberg, A. Hyva¨rinen, and F. Esposito, “Validating the indepen-dent components of neuroimaging time series via clustering and visual-ization,” Neuroimage, vol. 22, no. 3, pp. 1214–22, Jul. 2004.[87] N. N. Stone, L. Potters, B. J. Davis, J. P. Ciezki, M. J. Zelefsky,M. Roach, P. A. Fearn, M. W. Kattan, and R. G. Stock, “Customizeddose prescription for permanent prostate brachytherapy: insights froma multicenter analysis of dosimetry outcomes,” Int. J. Radiat. Oncol.Biol. Phys., vol. 69, no. 5, pp. 1472–7, Dec. 2007.107Bibliography[88] M. Keyes, D. Schellenberg, V. Moravan, M. McKenzie, A. Agranovich,T. Pickles, J. Wu, M. Liu, J. Bucci, and W. J. Morris, “Decline in uri-nary retention incidence in 805 patients after prostate brachytherapy:the effect of learning curve?” Int. J. Radiat. Oncol. Biol. Phys., vol. 64,no. 3, pp. 825–34, Mar. 2006.[89] S. Nouranian, M. Ramezani, I. Spadinger, W. J. Morris, S. E. Salcud-ean, and P. Abolmaesumi, “Learning-based Multi-label Segmentationof Transrectal Ultrasound Images for Prostate Brachytherapy,” IEEETrans. Med. Imaging, vol. 35, no. 3, pp. 921–932, March 2016.[90] I. Tosic and P. Frossard, “Dictionary Learning,” IEEE Signal Process.Mag., vol. 28, no. 2, pp. 27–38, Mar. 2011.[91] G. Sanroma, G. Wu, K. Thung, Y. Guo, and D. Shen, “Novel Multi-Atlas Segmentation by Matrix Completion,” in Int. Work. Mach. Learn.Med. Imaging, Sep. 2014.[92] M. Ramezani, K. Marble, H. Trang, I. S. Johnsrude, and P. Abolmae-sumi, “Joint sparse representation of brain activity patterns in multi-task fMRI data,” IEEE Trans. Med. Imaging, vol. 34, no. 1, pp. 2–12,Jan. 2015.[93] M. Aharon, M. Elad, and A. Bruckstein, “K -SVD : An Algorithm forDesigning Overcomplete Dictionaries for Sparse Representation,” IEEETrans. Signal Process., vol. 54, no. 11, pp. 4311–4322, 2006.[94] J. Tropp, “Greed is Good: Algorithmic Results for Sparse Approxima-tion,” IEEE Trans. Inf. Theory, vol. 50, no. 10, pp. 2231–2242, Oct.2004.[95] F. Sardanelli and G. D. Leo, Biostatistics for Radiologists: Planning,Performing, and Writing a Radiologic Study. Springer Science & Busi-ness Media, 2009.[96] M. Haber, H. X. Barnhart, J. Song, and J. Gruden, “Observer variabil-ity: a new approach in evaluating interobserver agreement,” J. DataSci., vol. 3, pp. 69–83, Jan. 2005.[97] S. Nouranian, M. Ramezani, I. Spadinger, W. J. Morris, S. E. Salcud-ean, and P. Abolmaesumi, “Automatic Prostate Brachytherapy Pre-planning Using Joint Sparse Analysis,” in Med. Image Comput. Com-put. Interv. MICCAI, vol. 9350. Springer International Publishing,2015, pp. 415–423.108Bibliography[98] G. E. Hinton, S. Osindero, and Y.-W. Teh, “A fast learning algorithmfor deep belief nets.” Neural Comput., vol. 18, no. 7, pp. 1527–54, Jul.2006.[99] L. Breiman, “Random Forests,” Mach. Learn., vol. 45, no. 1, pp. 5–32,Oct. 2001.109
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Information fusion for prostate brachytherapy planning
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Information fusion for prostate brachytherapy planning Nouranian, Saman 2016
pdf
Page Metadata
Item Metadata
Title | Information fusion for prostate brachytherapy planning |
Creator |
Nouranian, Saman |
Publisher | University of British Columbia |
Date Issued | 2016 |
Description | Low-dose-rate prostate brachytherapy is a minimally invasive treatment approach for localized prostate cancer. It takes place in one session by permanent implantation of several small radio-active seeds inside and adjacent to the prostate. The current procedure at the majority of institutions requires planning of seed locations prior to implantation from transrectal ultrasound (TRUS) images acquired weeks in advance. The planning is based on a set of contours representing the clinical target volume (CTV). Seeds are manually placed with respect to a planning target volume (PTV), which is an anisotropic dilation of the CTV, followed by dosimetry analysis. The main objective of the plan is to meet clinical guidelines in terms of recommended dosimetry by covering the entire PTV with the placement of seeds. The current planning process is manual, hence highly subjective, and can potentially contribute to the rate and type of treatment related morbidity. The goal of this thesis is to reduce subjectivity in prostate brachytherapy planning. To this end, we developed and evaluated several frameworks to automate various components of the current prostate brachytherapy planning process. This involved development of techniques with which target volume labels can be automatically delineated from TRUS images. A seed arrangement planning approach was developed by distributing seeds with respect to priors and optimizing the arrangement according to the clinical guidelines. The design of the proposed frameworks involved the introduction and assessment of data fusion techniques that aim to extract joint information in retrospective clinical plans, containing the TRUS volume, the CTV, the PTV and the seed arrangement. We evaluated the proposed techniques using data obtained in a cohort of 590 brachytherapy treatment cases from the Vancouver Cancer Centre, and compare the automation results with the clinical gold-standards and previously delivered plans. Our results demonstrate that data fusion techniques have the potential to enable automatic planning of prostate brachytherapy. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2016-06-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution 4.0 International |
IsShownAt | 10.14288/1.0305046 |
URI | http://hdl.handle.net/2429/58305 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2016-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2016_september_nouranian_saman.pdf [ 14.91MB ]
- Metadata
- JSON: 24-1.0305046.json
- JSON-LD: 24-1.0305046-ld.json
- RDF/XML (Pretty): 24-1.0305046-rdf.xml
- RDF/JSON: 24-1.0305046-rdf.json
- Turtle: 24-1.0305046-turtle.txt
- N-Triples: 24-1.0305046-rdf-ntriples.txt
- Original Record: 24-1.0305046-source.json
- Full Text
- 24-1.0305046-fulltext.txt
- Citation
- 24-1.0305046.ris
Full Text
Cite
Citation Scheme:
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>
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-0305046/manifest