使用灰階值與控制點加權的成本函數於不同病人間肺部非剛性對位

全文

(1)國立高雄大學資訊工程學系研究所. 碩士論文. 使用灰階值與控制點加權的成本函數於不同病人間肺部非剛性對位 Inter-Subject Human Lung Non-Rigid Registration with Weighted Intensity and Landmark Cost Function. 研究生：胡家祥撰指導教授：殷堂凱博士. 中華民國一百年七月.

(2) 使用灰階值與控制點加權的成本函數於不同病人間肺部非剛性對位指導教授：殷堂凱博士國立高雄大學資訊工程學系. 學生：胡家祥國立高雄大學資訊工程學系. 摘要. 在不同病人間肺部對位時，若僅做灰階值的對應，只能提供少量的資訊作為判斷依據。為了將轉換添加更多物理意義，氣管的分叉點被採用為控制點。因為其在不同病人間，雖然位置有不同，但是結構是相同的，如此便可替轉換增加實質意義。首先針對不同病人斷層掃描，找出胸部後將肺部取出，以適應性的閾值，搭配一個預先設定的信賴區間，達到快速且穩定的呼吸道擷取程序。並以擷取後的支氣管細化，以找出分支點。將不同病人間控制點配對後，進行影像對位。本論文提出一個加權的成本函數，使得影像對位得以同時對於兩個測量值收斂。成功得到肺部轉換，使影像可不僅將灰階值轉換對應，還可以拉近兩影像間肺部的內部結構。在五位病人交叉對位共二十五組配對，對照組僅能將平均控制點距離差靠近至 64.9 mm，本方法可拉近至 17.7 mm，改進 47.2 mm，改進幅度達到 73.8%。. 關鍵字：適應性區域成長、多解析度對位、加權成本函數、氣管分支點、肺部分區。. i.

(3) Inter-Subject Human Lung Non-Rigid Registration with Weighted Intensity and Landmark Cost Function Advisor: Dr. (Professor) Tang-Kai Yin Department of Computer Science and Information Engineering, National University of Kaohsiung. Student: Chia-Hsiang Hu Department of Computer Science and Information Engineering, National University of Kaohsiung. ABSTRACT Human lung image registration among patients provides little information if only intensity is mapped. Though the airway branch points are position varied, they are structure agreed. They can act as landmarks to register lungs. The thoracic section is extracted from CT, and the pulmonary region is segmented. An adaptive thresholding with a confidence interval is proposed to segment the airway. A weighted cost function is proposed to registering images with both similarity measure and landmark distance measure. It successfully registers lung images. 25 pairs within 5 patients are tested. The mean landmark distance of compared method is 64.9 mm. The proposed method achieved 17.7 mm that improved 47.2 mm, and it is 73.8% better.. Keywords: Adaptive region growing, Multi-resolution registration, Weighted cost function, Airway branch point, Airway segmentation. ii.

(4) Accknow wledg gement There are many peo ople I havee to thank for complleting this thesis. Thee most impoortant of all is my advissor Dr. Yin, Tan-Kai who w taught me m and direccted me. Un nder his superrvision, myy profession nal ability hhas grown, and a my attiitude and paassion reseaarching has innspired. I have to thhank to the thesis defennse committtee memberr, Dr. Chen,, Chia-Yen and a Dr. Hunaag, Wen-Chhen. They gaave me impportant advicces, and enrriched this w work. Of course, I would lik ke to thank tto my paren nts, brother and wife. T They supporrted me finishhing the graaduate progrram even thhere was a financial f issue. Then, I woould like to thank to m my colleagu ue of lab, Chen, C Cassiddy, and my y friend Hungg, Chi-Yanng. They gave g me uuseful ideass and moral supportss. Thanks to the comm munity of Innsight Segm mentation annd Registraation Toolkit, and Dr. M Marius Stariing and Dr. iir. Stefan Klein K of the elastix com mmunity, th hey gave me m lots of teechnical op pinions. Thannks to Chao, Joseph and Lee, I-Tinng, who gav ve me many y useful Engglish essay writing w tips. There are too t many people p I havve to thank k to. Please forgive mee that I cannot list everyyone. Thannks to all of you givinng me any supports in n the time I pursuing master degreee.. Chia-Hsian ng Hu National University U of Kaohsiu ung July, 100thh year of tthe Republiic Era. iii.

(5) Coonten nts 摘要要 ....................................................................................................................................... i Absttract................................................................................................................................. ii Acknnowledgement ............................................................................................................... iii Conttents ............................................................................................................................... iv List oof Figures...................................................................................................................... vi List oof Tables ..................................................................................................................... viii Chappter 1 Introdduction .......................................................................................................... 1 1.1 Prefacee .................................................................................................................... 1 1.2 Motivaation .............................................................................................................. 2 1.3 Purposee ................................................................................................................... 2 1.4 Compoosition ........................................................................................................... 3 Chappter 2 Backgground and Related Woorks .......................................................................... 4 2.1 Image Specificatio S on ................................................................................................ 4 2.2 Image Processing P .................................................................................................... 8 2.2.1 Mathematic M cal Morphollogy ......................................................................... 8 2 2.2.1.1 Basicc Operatorss ............................................................................... 9 2 2.2.1.2 Advaanced Operaators ...................................................................... 10 2.2.2 Image I Segm mentation ...................................................................................11 2.2.3 Image I Regiistration..................................................................................... 12 2.3 Relatedd Works ....................................................................................................... 13 2.3.1 Human H Airw way Segmeentation................................................................... 13 2.3.2 Lung L Regisstration ...................................................................................... 21 Chappter 3 Humaan Airway Segmentatio S on ........................................................................... 24 3.1 Thoraxx Extraction ................................................................................................. 26 3.2 Lung Segmentatio S n ............................................................................................... 27 3.2.1 Threshold T Selection S ................................................................................... 27 3.2.2 Connectivit C ty and Topoological Anaalysis ................................................... 28 3.3 Airwayy Segmentattion ........................................................................................... 30 3.3.1 Seed S Selecttion ........................................................................................... 30 3.3.2 Adaptive A Reegion Grow wing with Confidence Interval Thre resholding ...... 31 3.3.3 Leakage L Ch hecking ..................................................................................... 33 3.4 Postproocessing ...................................................................................................... 35 Chappter 4 Lung Registratio on............................................................................................... 36 4.1 Branchh Point Sets ................................................................................................. 37 4.1.1 Airway A Thin nning ........................................................................................ 37 4.1.2 Branch B Poin nts Identificcation ..................................................................... 38 iv.

(6) 4.1.3 Branch B Poin nt Sets Matcching ..................................................................... 38 4.2 Preregiistration....................................................................................................... 39 4.3 Registrration Architecture ...................................................................................... 40 4.3.1 Multi-Reso M lution........................................................................................ 41 4.3.2 Transform T .................................................................................................. 41 4.3.3 Optimizer O ................................................................................................... 42 4.3.4 Interpolator I r ................................................................................................ 42 4.3.5 Cost C Functiion............................................................................................. 43 Chappter 5 Experriment Resu ults and Disscussions ................................................................ 45 5.1 Airwayy Segmentattion ........................................................................................... 46 5.1.1 Thorax T Extrraction ...................................................................................... 46 5.1.2 Lung L Segm mentation .................................................................................... 47 5.1.3 Airway A Seg gmentation aand Postpro ocessing............................................... 47 5.2 Lung Registration R ................................................................................................. 50 5.2.1 Branch B Poin nt Sets ...................................................................................... 50 5 5.2.1.1 Airw way Thinninng ............................................................................ 52 5 5.2.1.2 Bran nch Points Iddentification and Matching ................................. 53 5.2.2 Registered R Images I ..................................................................................... 53 Chappter 6 Concllusions and Future Worrks ......................................................................... 61 6.1 Concluusions .......................................................................................................... 61 6.2 Future Works ........................................................................................................ 61 Reference ............................................................................................................................ 63 . v.

(7) List oof Fig guress Figurre 2.1 Streaak artifacts. (a) Stripes aacross the body b with 1.25 mm slicce thicknesss..... 6 Figurre 2.2 Noisee artifacts in n the bottom m of lungs. ................... . ........................................... 7 Figurre 2.3 A dem monstration n of dilationn and erosion n. Image is retrieved frrom Wikipeedia, the fr free encycloopedia, http://en.wikipeedia.org/wik ki/Mathemattical_morphhology. .......... 10 Figurre 2.4 A lettter B applieed with thinnning proced dure. Image is retrievedd from Wikiipedia, the free f encyclo opedia, http:://en.wikipeedia.org/wik ki/Topologiccal_skeleton. .11 Figurre 2.5 The concept c of image registtration. Imaage is retriev ved from thhe ITK manu ual. ............................................................................................................................................ 12 Figurre 2.6 A com mmon imag ge registratioon framewo ork. Image is i retrieved from the IT TK manuual. ............................................................................................................................... 13 Figurre 2.7 A hum man respiraatory system m. Image is retrieved r fro om Wikipeddia, the freee encyyclopedia, htttp://en.wik kipedia.org/w wiki/Respirratory_systeem..................................... 14 Figurre 2.8 Diam meter of airw way lumen iin the same patient. (a) Trachea, (bb) Main broncchus, (c) Innferior lobarr bronchus, (d) Anterom medial basall bronchus........................ 16 Figurre 2.9 A braain slice. ...................................................................................................... 22 Figurre 2.10 A thhorax slice. .................................................................................................. 23 Figurre 3.1 Flow w chart of thee proposed method.................................................................. 24 Figurre 3.2 Proceedure of thee proposed m method................................................................... 25 Figurre 3.3 Unwaanted cavitiies. ............................................................................................ 28 Figurre 3.4 Unwaanted cavitiies removedd. ............................................................................ 29 Figurre 3.5 The Suggest S area for retrievving seed po oints. ................................................... 31 Figurre 3.6 A Posssible brokeen airway w wall. ........................................................................ 33 Figurre 4.1 The Multi-resolu M ution registrration appro oach. Imagee is retrievedd from the ITK I manuual. ............................................................................................................................... 40 Figurre 4.2 The multi-resolu m ution registrration frameework. Imag ge is retrievved from thee ITK manuual. ............................................................................................................................... 41 Figurre 5.1 Extraacted thorax x section. (aa) 4 mm slicce thickness, (b) 1.25 m mm slice thickkness ............................................................................................................................. 46 Figurre 5.2 Segm mented pulm monary struccture. (a) 4 mm m slice th hickness, (b)) 1.25 mm slice s thickknes ............................................................................................................................... 47 Figurre 5.3 The leakage l into o lung. ....................................................................................... 48 Figurre 5.4 Segm mented airway in the prrogress of reegion growiing. (a) 4 m mm slice thickkness, (b) 1.25 mm slice thickness,, (c) 4 mm slice s thickness, (d) 1.255 mm slice thickkness, (e) 4 mm m slice th hickness, (f)) 1.25 mm slice s thickneess. ................................... 49 Figurre 5.5 Thinnned airwayss. (a) 4 mm slice thickn ness, (b) 1.2 25 mm slicee thickness. .... 51 Figurre 5.6 Thinnned airwayss with 4 mm m slice thick kness. (a) Co oronal projeection, (b) Axial A vi.

(8) projeection ............................................................................................................................ 51 Figurre 5.7 Thinnned airwayss with 1.25 mm slice th hickness. (a)) Coronal pr projection, (b b) Axiaal projectionn .................................................................................................................. 52 Figurre 5.8 Unwaanted airwaays branchinng. .......................................................................... 53 Figurre 5.9 Regisstered volum me in 3D. (aa) Source volume, (b) Preregistere P ed volume ..... 56 Figurre 5.10 Reggistered imaage in axial (left) and co oronal (righ ht) view. (a)), (b) Sourcee imagge .................................................................................................................................. 58 Figurre 5.11 Reaalizations try ying to findd the most su uitable weig ghting. (a) SSuccess ratio of differrent weighttings, (b) Im mprovementt of differen nt weighting g. ....................................... 59 . vii.

(9) List of Ta ables Tablle 2.1 CT sccans parameeters and speecificationss. ............................................................ 5 Tablle 2.2 Diameter of airw way lumen. ................................................................................. 17 Tablle 5.1 Iteratiions during region grow wing processs. ........................................................ 50 Tablle 5.2 Identiified branch h points withhin patientss. .......................................................... 54 Tablle 5.3 Comm mon branch points withhin patients............................................................. 55 Tablle 5.4 Compparison on laandmark di stance. ................................................................... 60 . viii.

(10) Ch hapteer 1 In ntrod ductio on 1.1 Prefacce In recent years, y radio ology is verry well dev veloped than nks to the improvemeent of technnology. Thee technique of image pprocessing has h also improved. Wiith the combined strenngth of thesse two techn niques, mo re details of o radiograp phy can be discovered d than ever before. Though it is common that X-raay images are presen nted in 2D,, it can alsso be preseented in 3D D ways usin ng Computeed Tomograaphy (CT) [15]. With the help of o CT, moree detail of spatial s relationship cann be provideed than trad ditional wayys. For exam mple, to rettrieve and reconstruct r 3D 3 anatomyy is easier with w the help p of CT. Therre are manyy applicatio ons that usses medicall image pro ocessing skkills, e.g. virtual v endooscopy [10]], nodule detection d [[48], organ n segmentattion and rregistration with anatoomy atlas [35, [ 41, 13]]. These skkills have utilized u man ny existing algorithmss, e.g. optim mal segmenntation [7], morphology m y operation [24] and im mage registraation [32, 45]. Thanks to the t promotiion of (Portaable) Electrronic Mediccal Records (P/EMR) by b the goveernment [1], medical images i cann be easily obtained with w P/EMR R by using g any portaable devicess on requestt, e.g. CDs, DVDs and d even thum mb drives etcc., in the form of Digittal Image annd Commun nication in Medicine (D DICOM) [4 43]. DICOM M is an indu ustrial standdard to form mat the disttribution annd transmisssion protoco ol of mediccal image so o that mediical imagess can be trransferred bbetween diffferent hosp pitals or m medical insttitutes withoout worryinng about the t inconsiistent file formats wiithin differrent scanneers or incom mpatible traansmission protocol p bettween scann ners and wo ork stations... 1.

(11) 1.2 Motivaation Intra-subject human lu ung registraation has become welll developedd and progrressed throuugh many difficulties in recentt years [50 0]. The sim milarity bettween each h CT exam minations is high, due to the same anatomical structure of patients soo that this task is not ttoo hard to accomplish h. The causee of the onlly difficulty y is subject to the positions, postuure and the breath hold ds of the subbject in diffeerent shooting time of C CT. It is muchh more chaallengeable to deal with w inter-su ubject imagge sets than n the intra-subject setts. Not only y the movem ments of thee bodies aree not the saame, but alsso the shapee of bodiess may be diifferent. In particular, the differen nt positionss, genders, races and aages may result r in som me tissue m might be thicker diameeter, and soome organs have higheer density. These vaariations m must be con ncerned in n the proceedure of image i regisstration. It is i simply no ot enough tto just conssider the sim milarity am mong the sub bjects for itt may mess up the expeeriment resuult by ignoring the variables.. 1.3 Purposse In order too overcome the lack oof characterristic differeences and tthe incompatible simillarity measuure, landmark distancce measure is introducced. By maaking use of o the innerr structure of the lung gs, the brannch points of o airway are a identifieed as landm marks. Thesse landmarkks are posittion varied but structu ure agreed among subbjects. Thuss, the Eucliidian distannce between n landmark sets could be considerred as a parrt of the meeasure of coost function. The combiination of similarity m measure and d landmark distance meeasure is a good soluttion and cann be used ass a cost funcction, becau use they aree both frequuently used as a the perfoormance inddex of a reg gistration prrocess. How wever, the balance of m making these two 2.

(12) meassurements as a a cost function f shhould be co onsidered carefully. c A weighted d cost functtion is propposed to strike a ballance betweeen similarrity measurre and land dmark distaance measurre. It attemp pts to mergee both advaantages of th hese measuurements wiith no muchh modificatiion.. 1.4 Compoosition This thesiss is organizeed as follow wed. First, in Chapter 2, backgrouund descrip ptions of raadiography and image processing are introdu uced. Relateed paper annd thesis aree also studiied and intrroduced. Th hey providee fundamen ntal techniq ques for furrther application. Thenn, in Chapteer 3, it desccribes the m method prop posed that reetrieves hum man airway y tree. In C Chapter 4, itt describes the main iidea of weiighted cost function ap approach off lung regisstration. Finnally, the exp periment reesult and disscussion, an nd the concllusion and future f workks are descriibed in Chaapter 5 and 6 respectiveely.. 3.

(13) Ch hapteer 2 Backg B groun nd and d Rellated Work W ks Before intrroducing th he contribut utions of th his thesis, th he backgroound and reelated workks are introdduced in thiis chapter. T These are reequired to be b known fo for understan nding the fo following chhapters. As mentionned in the previous chhapter, this chapter is organized to introducce the advaancement off radiograp phy and im mage processsing techniques respeectively. Reelated paperwork is inttroduced in section 2.3 .. 2.1 Image Specifiication In this worrk, there aree two groupps of 3D vo olumetric CT T scan sets . Each set of o CT data is a sequennce of 2D sllices. Indiviidual slices are stacked d into 3D sppace that forms a volum metric set. One of theese groups is low-dosee whole-body CT sets, which is provided by b the depaartment of Nuclear N Meedicine, Naational Chen ng Kung University U H Hospital, Taainan, Taiw wan, Republlic of Chin na (Taiwann). These data d sets were obtaineed by Siem mens's PET//CT scanneer. The averrage amount nt of slices in i each exaamination iss about 350 0 with the sslice thicknness equal to 4.0 miini-meter. Each E slice has 512 1.36. 512 pixeels of. 1.36 mm m spacing in width and lengtth. The slices are 122 bits preccision. monoochromatic images, bu ut stored in 116 bits resollution. It is acquired w with the follo owing param meters . 110 ~ 130 kVp an nd 50 ~ 1000 mAs;. . Spirall mode, 3 mm m collimatiion and 1.5 pitch; 4.

(14) . 4 ~ 5 mm m thickneess;. . 512. 512 matrix.. The other one o is low--dose wholee-lung CT scans, s which is retrieve ved from EL LCAP Publiic Lung Im mage Databaase, Vision and Imagee Analysis Group G / Intternational Early Lungg Cancer Action A Progrram (VIA/I--ELCAP) Public P Acceess Researchh Database. It is mainntained by the t I-ELCA AP Lab in tthe Weill Medical M Colllege of Coornell Univeersity, New w York, Uniited States of Americaa. The thorrax CT scan ns are obtaained in a single s breatth hold witth a 1.25 mm m slice thhickness wiith GE Lig ghtSpeed U Ultra. Their pixel spaciing is 0.76. 0.76 mm m, and storeed in 16 bitts grayscalee images. It is acquired d with. the fo following paarameters . 120 kV V, 80 mA;. . Helicaal mode;. . 1.25 mm m thicknesss;. . 512. 512 matrix.. The data seets provided d by the Nattional Chen ng Kung University Hosspital have. Tablee 2.1 CT scaans parametters and speecifications.. CT database d providers Specificatiion. ELCAP P Public Lu ung N NCKUH Imagge Databasee. Tubee voltage / current c. 110 1 to 130 kkVp / 50 to 100 mAs. 120 kVp / 80 mA. Scannning mode. Spiral S modee. Helical moode. Slicee thickness. 4 to 5 mm. 1.25 mm. Matrrix. 512 5. 512. 512 5. 5122.

(15) thickker slices thhan the oth her one. Thhe dataset acquired a by y the Corneell Universsity is relatiive thinner slices. Both h of them aare used in the t experim ments. The sspecification ns are summ marized in table t 2.1. Artifacts do d exist in both b groupss of scans, though t CT is a relativvely accuratte test [14]. They can be b roughly classified c ass 7 kinds. A.. Streakk artifact. Streaks cann be commo only seen arround materrials that blo ock most X X-rays, e.g. metal, m teethh and boness. It can be caused by undersamp pling, photo on starvationn, motion, beam hardeening, scattter, etc. Streak type off artifact offten occurs in the fosssa, or any metal implaants. It can be reduced by some addvanced recconstruction n techniquess [12]. B.. Partial volume efffect. It is some kind of blurring over sharpp edges. Scaanners may be not ablee to differentiate betw ween a smalll amount off high-densiity materiall and a larger amount oof lower deensity. Just like bones near cartilaage. The im mage processor of the scanner s triees to averag ge out. (a). (b). Figgure 2.1 Strreak artifactts. (a) Stripees across th he body with h 1.25 mm sslice thickness. (b) Teeth T block k x-ray resullt in streaks with 5 mm m slice thicknness. 6.

(16) the tw wo densitiees or structu ures. The deetail inform mation is losst, but the ccoarse lasts.. This can oonly be parrtially overccome by sccanners witth thinner slice thickneess. Two im mages with streak artiffacts are sho own in figurre 2.1. C.. Ring artifact a. It is probabbly the most commonn mechanicaal artifact. This T is the image of one o or manyy rings appeear within an a image. It is usually due d to a faullt of a detecctor. D.. Noise artifact. It is causedd by a low signal to nnoise ratio. This can commonly bbe seen in a thin slice thickness volume. It can also bbe seen wh hen insufficient powerr supplied to t the X-rayy emission tube to penetrate the obbject. It can n be seen in figure 2.2. E.. Motioon artifact. It occurs when w the object is shot w while movin ng. It results in blurringg and/or. Figure 2.2 Noise aartifacts in the t bottom of o lungs. 7.

(17) streaaking imagee. Motion blurring b cann be partiallly overcom me by Incom mpressible Flow Tomoography (IF FT) [23]. F.. Windm mill. When the detectors in ntersect the reconstruction plane, streaking aappearancess may occuur. With som me filtrationss or a reducction in pitch h can reduce this probllem. G.. Beam hardening. This can giive a cuppeed appearannce. This occcurs when there is moore attenuatiion in the ccenter of thee object thaan attenuatioon of the ed dge. It can be b easily sooftware corrrected by soome filters.. 2.2 Image Processing The imagee processin ng techniquues those are being used in this articlee are mathhematical morphology m operations,, image seg gmentation and a image rregistration n. The followings are thheir basic in nstructions of the generral usages and a functionns.. 2.2..1 Matheematicall Morphoology Mathematical morpho ology (MM)) is a techn nique and th heory. It is used to an nalyze and pprocess geoometrical sttructures. Itt is based on o Set theorry, Lattice ttheory, Topo ology and R Random fuunctions. Mathematical M l morphology can be applied nott only on digital d imagges, but alsoo on graphs,, surface meeshes, solid ds and otherr spatial struuctures. It can c be appliied on both continuous and discrette spaces wiith both bin nary and graayscale imag ges. Mathhematical morphology m was first prroposed to deal d with bin nary imagess in the 196 64 by Frencch, Georgess Matheron and Jean Seerra [24]. Most M of the MM M operatiions were developed in thee 1970s, e.g g. Hit-or-misss transform m, dilation, erosion, e opeening, closiing, 8.

(18) granuulometry, thhinning, skeeletonizationn, ultimate erosion, con nditional bissector, and otherrs. Mathematical morphology was thhen generaliized and exttended to grrayscale functtions and im mages in thee mid-1970ss to mid-198 80s. In 1986 6, it was furrther generaalized to a ttheoretical framework f based on coomplete latttices. This generalizatio g on brought flexibbility to thee theory whiich enablingg applicatio ons to a mucch larger am mount of strucctures, e.g. color c images, video, graaphs and meshes. m It forrmed the prresent MM usagees for image processing. Two of the basic operaations, dilattion and ero osion will bee introducedd in the following sectioons individu ually. Two oof the advan nced operations, closingg and thinniing will bbe introduced in the follow-up secctions.. 2.2.11.1 Basic Operatorrs There are two t basic operators o inn Mathemattical Morph hology. Theey are developed for bbinary imagges, but no ow they haave been ex xpanded an nd generalizzed to gray yscale imagges and com mplete latticees. They all reequire a kerrnel or struucturing eleement (SE). The SE iss a image with w a predeefined shape which pro obe an imagge. It is ofteen seen as a point or som me small im mages whicch is much smaller s than n the input im mage. One of the two basic operators iss dilation. It I uses the SE S to probee and expan nd the shapees containedd in the input image ass shown in equation e 2-1 1.. A⊕B. B⊕A. ⋃. ∈. B. (2-1). Erosion is the other on ne basic opeerator. It haas a similar idea to dilaation. It worrks as 9.

(19) equaation 2-2 shoowed. An example is showeed in figuree 2.3. It deemonstratess the usagee of dilation n and erosiion. This shhape is orig ginal in blu e. The dilated result is showed inn green, an nd the erodee shape is inn yellow.. A⊖B. ⋂. ∈. (2-2). 2.2.11.2 Advan nced Operators There are a lot of advaanced operaators in Matthematical Morphology M y. Two operrators have been adoptted in this work w are intrroduced bellow. Morphologgical closing g is a combbination of dilation d followed by errosion as sh howed in eqquation 2-3. The closing g operator iis used to reemove smalll holes in coontrast to. Figure 2.3 2 A demonnstration of dilation and d erosion. I Image is rettrieved from m Wikipediaa, the free en ncyclopediaa, http://en.wiikipedia.orgg/wiki/Math hematical_m morphologyy. 10.

(20) morpphological opening o whiich is used tto remove small s objectts.. A∙B. A⊕B ⊖B. (2-3). Thinning algorithm a or o morpholoogical skelleton is nott as naive as opening g and closiing operatorrs. It is an operator deesigned to extract e the medial linee of the shaape as show wed in figurre 2.4. It can be achieeved in sev veral ways. A hit-or-m miss transform is comm monly usedd to preservee the topolo gy of the sh hape.. 2.2..2 Imagee Segmen ntation Image segm mentation iss to divide a digital imaage into sub bsections. Th The goal of image i. Figure 2.4 4 A letter B applied witth thinning procedure. I Image is rettrieved from m Wikipediaa, the free en ncyclopediaa, http://en n.wikipedia..org/wiki/To opological_skeleton. 11.

(21) Figu ure 2.5 Thee concept off image regiistration. Im mage is retrieeved from tthe ITK man nual.. segm mentation is simplifying g or partitiooning the im mages, to maake the imagge much recoggnizable andd understan ndable. Precisely, it is a pro ocedure taggging pixels with con nsistent chaaracteristic. Like bounndary detecttion, it classsifies a pixxel that is belonging b to one objecct or anoth her by deteccting the color, brightness or textuure. There are lots of app plications uusing imagee segmentattion techniqque, e.g. tu umors locatting, tissue volume v meaasuring andd anatomy sttudy or diag gnosis usagees.. 2.2..3 Imagee Registrration Image registration is the t process that transfo orming an image i to annother coord dinate systeem. Input data can be images takeen from diffferent timee or differennt views. It must have some com mmon criteriia that bothh images haave, just lik ke the samee spot or ob bjects. The registrationn process triies to find a parameterr of the speecific spatiaal transform mation that bbest fit the target t imagee as shown is figure 2.5. For instaance, image registration n is to find a transform m parameter with transfformation T that transfo orms and m maps each piixel p in soource image to target piixel q in targget image.. 12.

(22) Figure 2.6 A commoon image reg gistration frramework. Imaage is retrieeved from th he ITK man nual.. The image registration n is usually treated as a framework k that has seeveral replaaceable com mponents as shown in fiigure 2.6. Itt has . a movving image which w is thee source imaage;. . a fixedd image wh hich is the taarget image;. . a metrric contains the cost funnction;. . an inteerpolator th hat interpolaates the pixeels not on th he grid;. . an opttimizer that provide thee strategy ho ow the metrric is minim mized;. . and a transformattion that maaps the mov ving image to t fixed imaage.. 2.3 Related Work ks Related papper and thessis those suppporting thiis work are introduced in the follo owing sectioons respectiively.. 2.3..1 Humaan Airwa ay Segmeentation n In the tradiitional way, the identiffying proceedure of airw way is mannually perfo ormed by a skilled raddiologist wh ho recognizzing local Regions R Of Interest (R ROIs) to fin nd out 13.

(23) Fig gure 2.7 A hhuman respiratory system. I Image is rettrieved from m Wikipediaa, the free en ncyclopediaa, http://en n.wikipediaa.org/wiki/R Respiratory_ _system.. wherre airway iss. The abnormal areas are then recognized th hrough the sscan series,, slice by sllice. Howevver, analyziing a large amount of series is a difficult annd huge wo ork. It takess lots of tim me, up to hou urs, to compplete the ideentification. The respiraatory system m is compossed with no ose, mouth, throat, airw way and lun ngs as figurre 2.7. It begins b at the t nose an and mouth, and termiinates at thhe alveoli. It is synonymously with w respirratory tract.. The air coming c from m external environmeent is inhalled from noose or moutth, and passses through h the pharyn nx into the ttrachea. Traachea 14.

(24) and bbronchi are the compon nents of airw way. The trrachea then separated innto left and d right broncchi at the caarina, which h is a cartilaaginous ridg ge at the low wer end of tthe trachea in the levell of the secoond thoracicc vertebra. T The main brronchi bran nch into bronnchioles, on ne for each lobe of eacch lung. Thee bronchiolees are subdiivided about more thann 20 times an nd go furthher into distal part of th he lobes. Affter division ns, the bron nchioles beccome shorteer and thinnner than beffore. They end up at tthe alveoli.. Alveoli arre anatomiccal structuree that form ms hollow caavities for gas exchangiing with thee blood [38]]. Human airrway is a trree-like tunnnel that air passes thro ough. It staarts from traachea and bbifurcates at a the carinaa which is rrecognized as the root node. Afterr the bifurcation, the trrachea is diivided into two t parts. T Those two are a called main m bronchii. Each bron nchus brancches into tw wo or more bronchioles b s, and each branch b poin nt is recogniized as a no ode in the trree. After branching b ab bout 17 tim mes or less, there are allveoli wheree gas exchaanges. Alveeoli are recoognized as th he leaf nodee of the treee [8]. It is an easyy task to seg gment and rrecognize th he trachea. The T tracheaa is located at the midddle of throatt and it is a rod shaped hollow tube as shown in figure 2..7 which co onnect from m mouth andd nose to maain left and right bronch hi. Althoough the seegmentation n of tracheaa can be deealt easily, the t segmenntation of aiirway becomes difficuult after seveeral bifurcattions. It is because b thatt the diametter of broncchiole becomes smalleer and smalller in the diistal bronchus. In figure 2.8, it shoows the diam meter of aiirway lumenn in differeent branchinng order. In n table 2.2, it shows thhe decremeenting diam meter of airw way lumen. And from tthe table, th he airway lu umen diameeter is obviiously decreeasing with each branching level. Moreover, with curren nt CT scannners, the qu uantity of slices s that llung occupiied is raiseed. And the resolution of images becomes higher h than ever beforee. The region of lungss across within w slicess increases from 80 to t a few hundreds. A According to the 15.

(25) (a). (b). (c). (d). F Figure 2.8 Diameter D off airway lum men in the same s patient nt. (a) Tracheea, (b) Main n bronchus, (c)) Inferior lo obar bronchuus, (d) Anteeromedial basal bronchhus.. advaancement off technolog gy, there is a significan nt differencce in how m many slicess that lung occupied. An automatic segmenntation way y is in neeed to deal with such large quanntity of slicees, and to replace this laabor-intensiive and timee-consuminng work. There is noo a fixed optimal o threeshold value, because the intensitty of the aiirway 16.

(26) Table 2.2 Diameter of o airway lum men.. Inferior Main M Trachea. Anterromedial. Pulmon nary. basall bronchus. alveoluss. lobar Bronchus B bronchuss Distal. Bifurrcation 0. 1. 2. 3 terminattion. Ordeer. 200 to 300 Diam meter. 1.98 cm. 1.78 cm. 1.59 cm. 0.87 cm c. m [46]]. lumeen is a little bit higher than t air inteensity. The range r of lum men intensiity over air is not a abssolute valuee. The intensity of air iss defined to o -1,000 Hou unsfield uniits (HU), an nd the intennsity of waater is defin ned to 0 H HU. Other intensity off tissue is aaccording to t the attennuation of X-ray. Thu us, as menntioned beffore, the optimal o thrreshold value is non-ppredeterminnable becau use of the arrtifacts menttioned in seection 2.1. Thus, theree are several ways to determine the thresho old. One off the metho ods is propoosed in ann interactive fashion [4]. Shorteening the operation o tiime withou ut big accurracy changeed, but it stiill requires an operatorr to deal witth the selecttion of cand didate voxeels. Nowadaays, basicallly, is to applly a manuallly selected fixed threshhold to acqu uire a roughh draft of thhe airway. The T selectedd threshold value is usually arounnd -1,000 to o -700 HU. Manually selecting one is a straightforw ward appro oach, becauuse the op ptimal thresshold differss from each people. Recent appproaches can n be classiffied into two o major typ pes: one is 33D approach h; the otherr is 2D and 3D hybrid approach. a Three-dimeensional typ pe is the appproaches th hose segmen nting the aiirway tree with w a 17.

(27) CT im mage seriess combined into a stackk of volumee. They are usually varried of 3D region grow wing algorithhms [34]. Many M of them m employ a classifier or manuallyy determinee seed pointts for a 3D seeded reg gion growinng procedurre [58]. This techniquee often resu ults in over--segmentedd, that is caalled the leeakage into lung paren nchyma, orr explosion. The airwaay wall is, sometimes,, possible bbeing brokeen in the low wer order bbronchi. It is i not only because off the limited resolutionn of the CT T scanner, but b also thee interferen nce of noisees and artifaacts such as partial voluume affect [20]. [ Hybrid typpe is the ap pproaches ccombining 2D 2 analysiss or reconsstruction an nd 3D segm mentation. It can be co ombined inn different orders o of th hese compoonents. Som me of them m perform 2D D filters to recover airrway walls, and then peerform 3D ssegmentatio on [5, 8, 100, 34, 59]. It can alsso perform in the rev verse orderr that meanns operating 3D segm mentation too identify th he main bronnchi and th hen 2D techn niques recoonstructing lower l orderr bronchi [33]. The segmeentation som metimes neeeds a predeffined thresh hold to deteermine whetther a voxeel belongs too airway or body whilee doing regiion growing g or segmenntation procedure. Becaause of the trachea, t it is not difficuult to be seg gmented. It is surroundded with mu uscle, fat, w water or othher tissues with higheer density th han air or greater thann -100 HU U. The trachhea is a ringg of airway wall filled w with air witth the intensity aroundd -1,000 HU U. The differrence betw ween the traachea and s urrounding tissue is significant. The intensiity of lungss is distribuuted from -1,000 to -6600 HU. Therefore, T itt seems nott too difficu ult to speciify a fixed threshold; however, h thhere is no a fixed value exists thaat suitable for f all casess. The airw way and thee lung havee similar inttensity, so that t they arre ambiguo ous to distinnguish. Although machine m leaarning technniques are useful u in so ome researcch areas witth the reduccing manuaal operationss, but the acccuracy of detecting d th he airway coonnectivity is not 18.

(28) goodd enough ussing a fuzzzy logic appproach [57]]. It segmen nts limited order of aiirway brancches. The most time-consu uming proccedure is to o find the best fittedd threshold.. The thresshold cannoot be retrieveed directly bby any prio or knowledg ge because iit may vary from persoon to persoon, and inffluenced byy artifacts mentioned in sectionss 2.1. The time consuumed by thhe processin ng of airwayy segmentattion depend ds on how thhe threshold d was choseen. An adaptivve way wass by propossed Kiraly et al. that tried t to solv lve the threeshold probllem [4]. Thhey tried to find the opptimal thresshold by ussing increm mental thresh holds, and they foundd most suitted thresholld value taakes place before the region gro owing proceess adaptedd to the en ntire lung. It is causeed by the leakage l of region gro owing proceedure that iddentifying outside o the aairway. The leakagge is also caalled as an eexplosion, because b the identified vvolume incrreased dram matically in a short perio od. Thoughh it actually catches the best one ass expected, it is a time consumingg algorithm just like ussing a sequeential search h algorithm m or a brute force one. In this workk, experimeents showedd that it perfformed as sllow as it seaarched. This algoriithm is not applicable while proccessing with h low resoluution scanss. The loweer the resoluution is, the fewer voxeels within each slice would be, andd also the higher h intennsity each voxels v are. More M speci fically, vox xels in the lung and airrway region n, not only have high intensity parts p such aas bones, vessels v and airway waall, but also o low intennsity parts liike the air in i the alveooli and airw way lumen. The T low inttensity parts will be hiigher than usual u while taking CT scanning and a sampling objects onn the fixed grids with low resoluttion due to aliasing a andd partial vollume effect,, and vice ve versa. The entire images willl looks mildd and smootth without detail d inform mation insteead of the aactual compposition. The averagee intensity becomes higher than normal or high 19.

(29) resollution ones. However,, the seed point withiin the trach hea has low w intensity. The intennsity in the bronchi b is much m higheer than the seed s point. The aliasingg will take place to sm mooth the entire scan n because there are vessels v and d micro veessels near each broncchus. There are not enough e griids to repreesent the im mage with low resolu ution. Moreeover, the distribution d of the vesssels and bro onchi is not uniform diistributed. Those T phennomena meentioned above, takes place in each e region n with diffferent influ uence. Therrefore, the thhresholds arre difficult tto discover while doing region groowing from m seed pointt to distal brronchi. Recently, Fabijanska F proposed a method that t takes the t usage oof the classified regioons with ann initial thrreshold [3].. It updatess the thresh hold whichh is going to be appliied by calcuulating the mean m valuee of the classified regio on. Having tthreshold arround the m mean value and plus an n interval. Iddeally, it con nverges while all airwaay candidatees are classsified. It trieed to choosse a better threshold quickly, q butt sometimess failed witth the exploosion of leaaking into th he lung entiirely. One reeason of thee failure of this algorithm is the m mean value of these vo oxels may nnot convergee due to thee incomplette airway wall w or the oover-estimaated step of choosingg a thresho old [33]. Moreover, M itt does not have termiinating connditions to terminate thhe procedurres beside th he convergeence of upd dating itselff. It often failed fa while dealing wiith the scan ns those hav ving thin sllice thickness. In otherr words, thee higher resolution thee scans are, the lower difference between th he air and nnon-air objeects would be. So the threshold is i easily over-updated exceed the most suitaable one, andd leads the segmentatio s on toward failure. fa In summarry, the meth hod of Kirally et al. pro oposed work ks fine withh high resollution scanss but workks slow with low re solution on nes, becausse of the sslow grow wth of thressholds beingg chosen. Itt increases 1 HU in eacch step. On the other hhand, the method of Faabijanska proposed p works well w with high resolution r scans, s but ooften fails while 20.

(30) dealiing low resoolution oness. The propoosed method d in chapteer 3 aims to combinee both advvantages off both methhods. It is a trade-off between sppeed, accurracy and su uccessfulneess. None of o the methhods mentiooned above perform onnly once reegion growiing with thee entire volume. Theyy are all inn an iteratiive fashionn. They alll require ex xtraordinaryy procedurres to accom mplish thee full segm mentation, e.g. anatomically an nalysis andd reconstruction. How wever, post-pprocessing is substituttable for eaach other allgorithms. H Hence, the main goal is to modify the regio on growingg part only, trying to ease e the com mputationall time and rrobustness.. 2.3..2 Lung Registra ation Image regiistration on n medical aapplication, especially y for brainn data, has been developed progrressively in n recent yearrs. Human brains b can be b registereed easily because of their rigid chaaracteristic. Though huuman brain is consists of soft tissu ues, e.g. greey matter, w white matteer and cerebbrospinal fluuid, luckily y, it is surrouunded by th he cranium, and is show wn in figure 2.9. The cranium iss solid tisssues those are hard and a difficullt to be beent. It not only consttraints the brain not to move eassily, but also preventss the brain being dam maged. Thuss, the shapes of brains are quite thhe same witthout much distorted, eexcept accid dental damaage or conggenital defecct. The bounndary of braain is circulaar and smoooth. Rigid regisstration show ws good ennough resultts registerin ng pairs of hhuman brain n and anim mal brain [17]. It is a compositioon of translations and rotations. C Combining rigid regisstration, scaaling and sh hearing, whhich is callled affine registration, r , shows a better b resullt [37]. Taking locaal deformation into connsideration, non-rigid registrations r s are propossed to 21.

(31) Figuree 2.9 A braiin slice.. deal with the deformable ability. Th The major non-rigid n registration r is the B-sspline regisstration which is a kern nel-based reegistration [11]. [ It put the t source aand target image i into iindividual finite f grids. It registers by moving each grid points. p There are still otheer kernel-bbased regisstrations liike Thin-PPlate-Spline and Elasttic-Body-Sppline registtration [39] . They aree proposed with advaanced usagees on deforrmed variattions. Lung is a soft and elastic organn filling witth air. It is much holloower than brain. b Therre have no many m boness beside ribbs surroundiing the lung g and show wn in figure 2.10. Withhout the connstraint of bones, b the vvolume and d surface haave huge diffferences am mong the luungs in everrybody. To get rid of the diffficulties off local defo ormations, researchers r at first atttempt applyying existedd methods of o brain-subj bjected regisstration to lu ung registraations. Rigid d way 22.

(32) Figure 2.10 A thorrax slice.. tried to overcom me this issuee but the ressult is still unable u to progress. Afteer that, non-rigid way shows a siggnificant im mprovement on lung reg gistration [56]. Making asssumptions and a constraiints may im mprove the performanc p ce of registration. Cao et al. proposed an approach tryiing to retain n the pulm monary voluume with a good resullt [27]. As mentionned in sectio on 1.2, a gaap does exisst in extendiing intra-subbjects operaations towaard inter-subbjects ones.. Not only the position n and statu us has to bee taken caree, but also the charactteristics amo ong people must be taaken into co onsiderationn. Ehrhardt et al. propoosed a non--linear land dmark-basedd approach trying to co onquer this issue with a fair resullt [21, 37].. 23.

(33) Ch hapteer 3 Huma H an Airrway Segme S entation The propossed method d, airway treee segmenttation, is inttroduced beelow [19]. There T are fo four stages of o this meth hod. The alggorithm is sh hown in the following fflowchart, figure f 3.1. JJust like anny other airw way segmenntation algo orithms, the first two sttages and th he last stagee are quite thhe same. Th he third stagge is the corre part of th his algorithm m.. Figuree 3.1 Flow cchart of the proposed method. m 24.

(34) Inpu ut:. CT C scan. Outp put:. Mask M of airw way. Steps: 1.. thoorax = ExtraactThorax ( volume ). 2.. lunngs = SegmentLungs ( thorax ). 3.. traachea = Seg gmentTracheea ( thorax ). 4.. seeeds’ = seedss = SelectSeeeds ( trach hea ). 5.. wh hile ( ! leakaage ( lungs’’ ) ). 6.. thresholld = averagee ( seeds’ ). 7.. lungs’ = threshold ( lungs ). 8.. for each h voxel in luungs’. 9.. if connected ( seeds’ , vo oxel ). 10.. seeds’ = add ( seed ds’ + voxel ). 11.. maask = region nGrowing ( lungs , seed ds , thresholld ). 12. maask = postPrrocessing ( m mask ) Figure 3.2 Proceddure of the proposed method. m. As shown in i figure 3.2 2. The entiree algorithm m is describeed and discuussed step by wing. step as the follow . Read the t entire volume slicee by slice. . Extracct the thorax x regions. . Segmeent the lung gs. . Segmeentation of trachea andd bronchi . S Select seeds. . C Calculate the average oof seed voxeels. . 3 adaptive region grow 3D wing with mean m and co onfidence innterval. . U Update the mean m valuee with the average a of segmented vvoxels whicch are b belongs to airway. . C Check leakaage into lungg 25.

(35) . If not leeaking, go bback to regiion growing g step. . If leakiing into lunng, return the t result of o region gr growing witth the previou us mean vallue. . Post processing p. . Outpuut result. The steps inn each operration are deescribed in the followin ng sections... 3.1 Thoraxx Extra action First of all,, in order to o achieve thhe goal of seegmenting the t airway ttree, the loccation of thhe airway must m be know wn. The airrway is up to noses an nd down to diaphragm as in figurre 2.7. The interesting part is loc ated within n the lungs. The lungs are parts of o the chestt. The regioon of intereest is locate d at the cheest. It is neecessary to extract the chest sectioon from a whole-body w y scan, but iit can be skip and no need n if a whhole-lung sccan in somee databases is used. Thorax exttraction can n be done m manually with ease. Th he thorax seection is ussually locatted from abbout 100th slice s to 1900th in an ordinary 4 mm m upper-boody CT scan n. An extraaction contaains from sh houlders to diaphragm in axial ord der. Shouldders and necck are obvioous distinguuishable. Th he bones o f both armss show up after shouldders. The end of the luungs is a liittle bit morre difficult but still disstinguishablle by eyes. The intensity in the luungs is low w because lu ungs are fillled with airr, and the in ntensity of tthe organs under u the luungs is highh because liiver, stomacch and intesstines are co omposed wiith fat, tissu ue and water those intennsity are deefinitely greeater than air. There is allso a fuzzy logic approoach that ussing body width, w volum me of cavity y and heighht of candidate area as a criteria too extract th he thorax section [2]. Though it is an autom matic algorrithm, it is useless in the situatio on of alread dy extractedd scans and d less 26.

(36) accurrate than manually extrracted one. In this worrk, thorax section s is m manually ex xtracted in order to foccus on the main proceedure of airrway segmentation.. 3.2 Lung Segmen S ntation After chestt extraction, the secondd stage is to o segment respiratory r ssystem in th horax sectioon that incluudes lungs and airway for further segmentatio on. There are two t fully feeatured algoorithms thatt already well developeed. One of these two aapproaches,, Hu et al. proposed p ann optimal seegmentation n using optim mal thresho olding [49],. another. one. proposed. bby. Sluimeer. et. al... which. performin ng. a. segm mentation-byy-registratio on scheme [[18]. Both of o them are excellent w works to seg gment humaan lungs, annd they are both b accepttable for dirrect adaption n in this theesis. The methood Hu et al.. proposed has lower complexity, because m machine leaarning mechhanism and image regisstration havve not been taken used [6, 37]. It iss so fast to apply somee simple thrresholds insttead of largge amount of decision rules being aapplied. The instrucctions of thiis algorithm m are introd duced below w. After all, the later part of origiinal algorithhm is deprecated and no need to t apply, th hat eliminatting tracheaa and broncchi, becausse those arre what inn need in this thesis. In this thhesis, cascading operaations are going g to app ply from herre.. 3.2..1 Thresh hold Selection Equation 3-1 shows th he optimal tthreshold. Itt is an iterative proceddure selectin ng the thresshold T. Thhe μ and μ are thee averaged intensity of o the bodyy and non--body voxeels respectivvely. 27.

(37) Figure 3..3 Unwanted cavities.. T. (3-1). It is assumed that the image contaains only tw wo types of voxels, boddy and non--body voxeels. Lungs and other low densityy tissues classified c ass non-bodyy voxels. Voxels V withiin the body those have high densitty are classiified as body voxels. The initial threshold is i set to -1,0000 HU because the in ntensity of air is predeefined as –11,000 HU. The T new threshold in tthe next iteration is thee average oof mean inteensity of booth body andd non-body y voxels withhin the chesst.. 3.2..2 Conneectivity and a Topoological Analysiss There exist some smaall componnents those are not con nnected to the major parts, lungss. Three dim mensional connected c ccomponent labeling is applied to iidentify thee lung 28.

(38) Fig gure 3.4 Unnwanted cav vities remov ved.. regioon. Componnents with the two laargest volum mes are thee left and rright lungs. The regioon, with tinny volume and is noot connecteed to the lungs, will be discard d and elimiinated. Interior caavities in th he lungs arre unwanteed. Vessels filling witth blood, which w consiist of wateer and blood cell withh high denssity, are rellatively higgh intensity y than surroounding air. These will be removedd during thee thresholding procedur ure. This phenoomenon is sh hown in figgure 3.3, and d can be eassily recoverred by finding its comppliments [222]. At first, this algoritthm, find th he complimeent region oof the segm mented lungss. Next, it finds the complimennt region which w is connected c w with the laargest compponent. Thee operation can c be simpplified as eq quation 3-2 and figure 33.4.. R ← con nncomp R , 1 29. (3-2).

(39) 3.3 Airwayy Segmentation n Human airrway segm mentation iss implemen nted using three-dim mensional seeeded regioon growing with a seg gmented thhorax CT sccan. It appllies region growing with w a speciific criterioon from a seed point w which is loccated at thee center of the tracheaa, and then applies a ruun time deteermined threeshold insteead of a fixeed one. The detaileed individu ual proceduures of segmentation are describbed in follo owing sectioons respectiively.. 3.3..1 Seed Selection S n The seed point p is stro ongly suggeested to be selected fro om the centter of the aiirway lumeen. It is beccause that th he intensityy decreases from airwaay wall to tthe center of o the airwaay lumen. The T airway wall w has higgh intensity y relative to the inside oof airway lu umen. Thosse tissues haave significaantly higherr density thaan air. It is trivial to find a circular annd low intensity area in tracheaa locating at a the superrior slice rather r than finding broonchi in th he remainin ng slices. Itt is becausse the trachhea has largeer diameter than bronchhus. It is importtant to selecct a seed poiint near the center of aiirway lumen en. The closer the seed point to thee lumen cen nter is, the lower the in ntensity retrieves. If a seed is obttained near the lumen wall, it wiill have higgher intensiity than thee center one ne, due to partial p volum me affect. It I is highly risky resultting in leak king into lun ng while peerforming region grow wing proceddure. Voxelss belongingg to the walll are easily y miss-classsified becau use of the ssimilar intennsity retriev ved. The traachea is surrrounded by y body voxeels in contrrast to the bbronchus suurrounded by y air or ves sel. It rarely y leaks whiile processinng the trach hea. It usuallly leaks tilll processing g to the bronnchi which are surroun nded by alveeolus. 30.

(40) Figure 3.5 5 The Suggeest area for retrieving seed s points.. While dealling with thin slice scan ans, the salt and pepperr noise is m much easier being encountered relaative to the thick ones. It is possib ble retrieving g a seed poiint that is a noise pointt with muchh higher vallue than its own neighborhoods. The T solutionn is smooth h with a meedian filter instead i of mean m filter, bbecause thee mean valu ue can be eaasily raised, even by a single whiite noise. Another A plaausible way y is better to t add morre than one seed pointts. Adding more m than one o seed pooints will reesult in redu ucing the innfluence of noise pointts. Seed poiints are set in i the positiion of tracheea shown in n the figure 3.5.. 3.3..2 Adapttive Regiion Grow wing witth Confid dence In nterval Thresh holding The main idea i of thiss approach is segmentaation by reg gion growinng [61]. It starts with an initial thhreshold that is the inttensity of th he seed poin nt. Its intennsity is as lo ow as 31.

(41) air. T The lower inntensity, as low l to -1,0000 HU, willl lead to bettter segmenttation resultt. The threshold value going g to appply for regio on growing is determinned by the mean valuee of those voxels v alreaady being cclassified ass part of airrway tree aand a confid dence intervval. The mean m value is multiplieed with a manually m asssigned connfidence intterval whicch is assigneed by being g experimennted. While doing regio on growingg, more and more voxeels are classiified as partt of airway ttree. After addinng any vox xels join thee airway treee, it will result r in uppdating the mean valuee those voxeels classifieed as airwayy. By doing so, s the threeshold valuee becomes higher and d higher. Thhe growth of o the increement of meean value will w slows ddown in a feew iteration ns, because of the addittional toleraance of threeshold valuee. The meann value is multiplied m with w a confiddence interv val. It is hoow this systeem slows itss growth annd converges. There existts a fatal flaw that theere is no en nough practiical terminaation ending g this proceedure. It terrminates by only waitinng for itselff convergence. It needss to predefin ne the amouunt of loopss that is goiing to iteratte. It breakss the loop iff the predefi fined threshold is reachhed, and resstricts the iteeration. Preeventing thee loop from over-iterateed. In the firstt time of region r grow wing iteratio on, the mean value oof seed poin nts is evaluuated as shoown in equaation 3-3.. (3-3). The threshhold value is i evaluatedd by the mean m value of classifieed region with w a confi fidence interrval as show wn in equatiion 3-4.. 1 32. 1% 1. (3-4).

(42) Figu ure 3.6 A Poossible brok ken airway wall.. It has beenn tested thaat one perceentage of co onfidence in nterval is thhe most suitable one. The threshhold value used u for reggion growin ng is the meean value m multiplied by y one minuus confidencce interval, which is 999%. While iteraating, the mean m value of classifieed region is i recalculat ated as show wn in equaation 3-5. Thhe threshold d is then uppdated. It means m the th hreshold havve 1% tolerrance. Larger interval will resultt in the le akage, and d the smaller one willl result in poor perfoormance.. (3-5). 3.3..3 Leakaage Checcking Due to the limited resolution, disstal bronchi wall are so ometimes brroken after some bifurrcations. Onnce it breaaks, the reggion growin ng procedu ure previouusly describ bed is highlly possible classifying g those voxxels with low w intensity in the alveeoli or bron nchus otherr than the processing one. o It is so called leakage or explosion into llung. As a result, r 33.

(43) it is difficult too identify th hose lower order bron nchi. A posssible brokenn airway wall w is show wed in figuree 3.6. There are two t ways to o deal with tthis issue. One O is to sm mooth or to downsample the entiree image. It makes the airway a walll linked, bu ut loses preccession. It iss also possib ble to fill uup the entirre airway lu umen, makking the reg gion growin ng procedurre impossib ble to reachh lower order bronchi. Another waay is to low wer the thresshold. With lower thresshold, theree will be feewer chancces to leak into lung, but still no chance to identify distal airwaay. Therefore, the region growing prrocedure has to make sure s that theere is no leaakage occuurred, and prrevents the threshold bbecoming to oo high. Fallling back th the threshold and gettinng the best one before leakage. The leakagge can be deetected by m monitoring the volumee of classifi fied region. Once the leeakage occuurred, the volume increeased dramaatically. Thee intensity oof voxels ou utside the aairway is disstributed in a small rannge, and therrefore it is easy e to expllode to the entire e lung.. In this woork, leakage detectionn is implem mented by monitoringg the segm mented regioon with a sppecific thresshold volum me. It tries to o catch the growth g of cclassified vo olume and qquantities how h dramatiically it gro wn. If an ex xponential volume v groowth of classified regioon is detectted as equaation 3-6, a leakage iss identified,, and then stops the region r grow wing processs.. (3-6). 34.

(44) 3.4 Postprrocessin ng Any tissuee in the natture world is rarely to o be seen with w sharp oor zigzag edges, e bounndaries andd surfaces. Airway iss an organ with smo ooth surfacee likewise.. The segm mented airw way that produced p byy previouss procedurees may haave cornerss and fragm ments. In order to t reconstru uct the bouundary of airway waall and fulffill the cav vities, morpphology cloosing operator is applieed. After app plying this operator, it would get better b visuaal effect andd reasonablee segmentattion.. 35.

(45) Ch hapteer 4 Lung L Regisstratiion In order to overcom me the laack of inteer-subjects human luung registration, splinne-based traansformation n is used too solve the inconsistent i t shape infoormation; aiirway brancch points arre used to matching m thee structure of o lungs. Though thee similarity measure is used by default, it is not n sufficiennt to registerr with only this measurrement. It ignores the pphysical meeaning of im mage presennted. It conssiders only the mappinng of gray value, v and caausing the registration r become meeaningless. The airwayy branch poiints are adoopted as land dmarks, beccause the aiirway determ mines how the shape of lung will w be. Thee airway is the inner structure oof lungs. It exist regullarity of thee airway bifu furcates. Theerefore, airw way branch points havee to be iden ntified as lanndmarks. In order to take the ad dvantage of the local strructure, it has h to take uuse of the branch pointts. Formingg the landm marks as a measuremeent in behaalf of the qquantificatio on of brancch points iss the way to o take brannch points as a a cost fu unction. Thuus, the Eucllidian distaance withinn landmark sets couldd be consid dered as a measurem ment of thee cost functtion which is i going to be b minimizeed. In this worrk, it adoptss the methodd Li et al. proposed p which tries too register human lungss by Thin-P Plate-Splinee transform mation [6, 28]. 2 It miniimizes all m measurement by exhaaust search,, and stopss at a pre defined am mount of itterations. IIt is a com mmon impleementation of image reegistration, but loses effficient. It is an impportant issue of the com mbination of o both meaasurements. A weighted d cost functtion is propposed for th he sake of kkeeping botth advantag ges of each measure. In n this workk, the best weighting of similarrity measurre and land dmark distaance measu ure is introduced. It trries to find the t best suiit weighting g by trial an nd error, andd trying to find f a 36.

(46) robust weightingg with good d performannce.. 4.1 Branch h Point Sets To take usee of the braanch point sset which reepresent the inner struccture of the lung, the aairway brancching must be identifieed. In the prevvious work in chapter 3, the airway is successfully segm mented. It is i not readyy for a direect usage. Because B thee airway is a tree-like air tunnels,, the decision of seleccting a voxeel represent as the brancch point is an a issue neeeds to overccome. The proceddures of colllecting brannch points are a describeed in the foollowing secctions in deetail.. 4.1..1 Airwaay Thinn ning The airwayy segmented d, is a tree-llike structurre. There arre airway w wall and lum men in each tunnel. Thee cross sectiion of the aiirway is a flat f lumen, so s that it is ddifficult to find f a singlle point reprresent the branch pointt [31]. In order too retrieve the mediall line of th he airway, a thinningg / skeleton nizing algorrithm is in need n [30]. Itt helps the ooperator gettting the ind dex of the brranch point. The thinninng procedu ure is to prooduce a skeleton of th he airway. It will be a one voxeel wide struccture, to easse the difficculty of retrrieving a sin ngle point th that represen nt the brancch point. The reasonn not to app ply pruning algorithms after thinniing proceduure is to sav ve the proceessing timee. A postprrocessing likke morphology closin ng or smoooth filters before b thinnning proceddure may better b help than prunin ng, becausee the pruniing result is i not recovverable. If a branch is miss-pruned m d, a branch point will be b lost. 37.

(47) In this worrk, the method Homannn et al. prop posed is ado opted. Homaann et al. already contrributed this implementation to thee Insight Jou urnal which h is an open access and d open peer--review on-line publicaation in meddical image processing and visualiization [16].. 4.1..2 Brancch Pointss Identifi fication Branch points are eaasy to identtify manuallly. After the thinningg proceduree, the imagge now is a skeleton of o the airw way. All thee identification proceduure has to do is traveerse the airrway tree. It I can be B Breadth-firsst search (B BFS) or Deepth-first search s (DFS S). Mori et al. proposed an a automaticc branch po oint retrievin ng algorithm m [29]. But there is onne thing haas to consid der, the airw way might be miss-seg gmented. Itt causes loo oping issuee and unwannted false branching b [660], so that Tschirren et e al. tried tto overcom me this issuee [25]. In this worrk, the brancch points arre identified d manually on purposee of accuracy and stabiility.. 4.1..3 Brancch Point Sets Maatching There are lots of bifurcations in thhe airway. In I order to find f commoon branch points, only the major branch b points are consiidered. Twenty-thrree branch points p with low branch hing order are a selected for recognizing. Theyy are locateed at trachea or carina , right main n bronchus, 1st divisioon of right main broncchus, right bronchus 1, right bronnchus 3, righ ht intermed diate, right bbronchus 6, right broncchus 4 + 5, right broncchus 7 + 8, right broncchus 9 + 10, right bronnchus 9 and d right broncchus 10. Onn the other side, they are located d at left maiin bronchuss, 1st divisiion of 38.

(48) left m main bronchhus, left upp per divisionn, left lingu ular division n, left inferiior, left bron nchus 6, lefft basal, leftt bronchus 8, 8 left broncchus 9 + 10,, left bronch hus 9 and leeft bronchuss 10. In this worrk, the airwaay branch p oint sets aree matched manually, m coonsidering errors e may be encounntered in prrevious proocedures. Th hough auto omatic algoorithms do exist, theree are still roooms for refiine as menttioned in preevious sectiion.. 4.2 Prereggistratio on In order to t speed up p the entirre process, a preregisstration is adopted. Image I regisstration is a time consu uming imagge processiing techniqu ue. If imagges have rou ughly alignned, preregiistration reeduces timee spent by the later registrationn which is high compputational expense. e In this worrk, a rigid-b body transfoormation is employed. It consists oof only rotaations and ttranslations. It is one of o the simpllest transforrmations tho ose do not rrequire too much compputational resources. There are six param meters to deescribe the rigid-body y transform mation matrrix as equaation 4-1.. M. TR. (4-1). The rigid-bbody transfo orm matrix iis equal to combining c translation t ttransform matrix m and rrotation trannsform matrrix, where. 1 0 T  0  0. 0 1 0 0 39. 0 q1  0 q2  1 q3   0 1. (4-2).

(49) and. 0 0 1  0 cos q  4  sin  q4  R  0  sin  q4  cos  q4   0 0 0. 0   coss  q5   0   0 0    sinn  q5   1  0. 0 sin  q5  0   cos  q6  sin  q6   1 0 0    sin  q6  cos  q6  0 0 cos  q5  0   0  0 0 0 0 1 . 0 0  0 0 1 0  0 1. (4-3). 4.3 Registrration Architec A cture The propossed method is designedd in order to o compete with w the metthod Staring g et al. propoosed [42], the components, multti-resolution n pyramid, optimizer and interpo olator remaain the sam me as consttraints and fixed facto ors. The co ost functionns are contrrolled factoors those arre different from the ccompeting method. m On ne of the traansformatio ons is added because itt is required d by the corrresponding cost function.. Figure 4.1 1 The Multii-resolution registration n approach. Imaage is retrieeved from th he ITK man nual.. 40.

(50) Figure 4.2 The multi-rresolution registration r frameworkk. Imaage is retrieeved from th he ITK man nual.. 4.3..1 Multi--Resoluttion In this worrk, a multi-rresolution aapproach is adopted. Itt is a coarsee to fine method that rregistering from f low reesolution to original ressolution forr the sake off speeding up u the entiree registratioon as shown n in figure 44.1. A Gaussiann blur filterr is used to build the multi-resolu m ution pyram mid. The oriiginal imagge is blurredd then subsampled to a lower level pyramid.. Repeatingg this step until u a user specific pyyramid levell reached. A multi-reso olution fram mework is sshowed in figure f 4.2.. 4.3..2 Transfform Betke et all. proposed a lung rigiid registratiion using Itterative Cloosest Point (ICP) way [37]. Now wadays, rig gid registrattion is nott good eno ough for luung registration. Non--rigid transfformation iss in need. According to the meethod Rex Cheung ett al. propossed, TPS iis just fitteed for landm mark-basedd registratio on [40]. T PS transfo ormation do oes not reqquire too many param meters. Onlly source an nd target lanndmarks arre required. Therefore, there is no need 41.

(51) to finne tune any parameterss. In this worrk, the TPS implementtation is bassed on the method m Davvis et al., Brooks and A Arbel propoosed [39, 47 7]. The transsformation is i illustrated d in equatioon 4-4.. ∑. (4-4). A postregiistration is required to fitting the shape informatioon by B-S Spline transsformation. It is because that the TPS only tries t to match the landdmarks. B-S Spline regisstration putss both imagees on finite grid which is denser th han landmarrks. As the result, it is nnecessary too adopt B-S Spline registtration to match m the bo oundary. B-SSpline is sh howed in eqquation 4-5.. ∑. ∈. (4-5). 4.3..3 Optim mizer In this work, Adaptivee Stochasticc Gradient Descent D (AS SGD) optim mizer is used d [51]. It is a stochasttic gradientt descent ooptimization n method for f image rregistration with adaptive step siize predictiion. In conntrast to Reegular Step Gradient D Decent (RS SGD), AGS SD has an addaptive step p size whichh helps the convergent c speed fasterr than ever.. 4.3..4 Interp polator Trilinear innterpolation n is used inn this work.. It approxiimates the vvalue withiin the locall axial data on o lattice po oints. It is a faast algorith hm with loow computtational complexity aagainst B-sspline 42.

(52) interppolation annd windowed sinc intterpolation, and resultts enough good qualiity of interppolated voxxels. It is a little bit sllower than Nearest Neeighbor (NN N) Interpolation, but pproduces bettter result.. 4.3..5 Cost Function F n The cost fuunction is th he most im mportant and d the core part p of this tthesis. The main contrribution is to t provide a practical weighting for f combiniing both sim milarity meeasure and llandmark diistance meaasure. The similaarity measurre is shownn in equatio on 4-6. It is the gray value diffeerence withiin images. Mean of Squared S Diifferences (MSD) is employed foor the simiilarity meassure in this work. w The MSD M is defi fined as,. MSD T ;. ,. Ω. ∑. (4-6). ∈. wherre Ω is the lung region n and μ is tthe instancee of transforrmation. mark distancce is showeed in equattion 4-7. Itt is aimed tto minimizze the The landm distaance of twoo point sets with knoown corresp pondence. It makes tthe local aiirway struccture matcheed. It compu utes the disttance from target landm marks to traansformed source landm marks.. LDM M T;. ,. ∑. x. T x. (4-7). The goal of o this thessis is to meerge two measuremen m nts as a cosst function by a weigghting. A weeighted costt function iss shown in equation e 4-8 8, 43.

(53) C T ;I ,I. ∑. ∑. ω c T ;I ,I. wherre C is thee weighted cost functiion, and respeectively. Each. and. is arrranged withh a weightin ng. (4-8). aree target andd source im mages for to otal N equal al to two.. The cost fuunction can be expandeed as equatio on 4-9.. C T ;I ,I. ω. ω. 1 Ω ∑. ∈. x. T x. (4-9). While norm malizing th he measureements, the maximum m and minim mum should be discoovered. Thee maximum m of MSD is equal to o square off twice maxximum inteensity availlable in DIC COM, and th he minimum m is zero. After A perform ming imagee registration, the imagges are matcched. Thereefore, the sim milarity is high. h If the MSD is noormalized, it i will be diistributed neear zero. Th he LDM is ssimilar to MSD. M In this casee, cost functions are w weighted witthout normaalization. Thhe normalizzation coeffficient is im mplicit and taaken place bby weightin ng directly.. 44.