一種基于形狀約束的膝關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取方法
【專利摘要】一種基于形狀約束的膝關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取方法，本發(fā)明涉及基于形狀約束的膝關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取方法。本發(fā)明為了解決現(xiàn)有半月板分割方法計(jì)算效率與處理效果的不平衡、利用單一灰度信息使處理效果不理想、需要人工設(shè)置初始點(diǎn)、難以實(shí)現(xiàn)自動(dòng)圖像序列批處理、提取效果受成像條件影響等問(wèn)題。具體過(guò)程為：一、ROI定位，選定初始片層組；二、進(jìn)行優(yōu)化，得到優(yōu)化后片層集合；三、對(duì)優(yōu)化后片層集合中的片層的ROI進(jìn)行閾值自動(dòng)選擇與半月板目標(biāo)提取；四、根據(jù)提取出的半月板目標(biāo)在優(yōu)化后片層集合中進(jìn)行半月板原始信息提取，得到包含分割出的半月板原始信息的ROI的片層。本發(fā)明應(yīng)用于半月板提取領(lǐng)域。
【專利說(shuō)明】
-種基于形狀約束的膝關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取 方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明設(shè)及基于形狀約束的膝關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取方法。
【背景技術(shù)】
[0002] 現(xiàn)有的膝關(guān)節(jié)半月板提取方法有基于目標(biāo)邊界信息的的分割方法，包括邊緣檢測(cè) 算子、區(qū)域生長(zhǎng)、分水嶺算法等。此類方法施行起來(lái)較為簡(jiǎn)單、運(yùn)算量小，但利用的信息較單 一、無(wú)法結(jié)合醫(yī)學(xué)經(jīng)驗(yàn)知識(shí)，目標(biāo)分割效果差;有基于輪廓模型的方法，包括snake模型及其 改進(jìn)算法，如主動(dòng)輪廓模型、測(cè)地線主動(dòng)輪廓模型等。此類方法處理效果較好，但計(jì)算量大、 模型演化速度慢、初始條件需要人工設(shè)置且對(duì)最終處理結(jié)果影響大，同時(shí)難W實(shí)現(xiàn)圖像序 列批處理;有基于灰度信息的提取方法，包括紋理目標(biāo)分割、闊值處理等。此類方法應(yīng)用了 更多的圖像特征及醫(yī)學(xué)經(jīng)驗(yàn)知識(shí)，可W較好的平衡計(jì)算效率與處理效果，但存在受成像條 件影響產(chǎn)生誤提取的現(xiàn)象，如ROI定位受成像位置影響、關(guān)節(jié)軟骨被誤提取為半月板目標(biāo)等 現(xiàn)象。

【發(fā)明內(nèi)容】

[0003] 本發(fā)明為了解決現(xiàn)有半月板分割方法中存在的W下問(wèn)題:計(jì)算效率與處理效果的 不平衡、未結(jié)合醫(yī)學(xué)經(jīng)驗(yàn)知識(shí)利用單一灰度信息使處理效果不理想、需要人工設(shè)置初始點(diǎn)、 難W實(shí)現(xiàn)自動(dòng)圖像序列批處理、提取效果受成像條件影響等問(wèn)題，從而提出一種基于形狀 約束的膝關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取方法。
[0004] 上述的發(fā)明目的是通過(guò)W下技術(shù)方案實(shí)現(xiàn)的：
[0005] 步驟一、ROI (感興趣區(qū)域)定位，選定初始片層組Sselected ;
[0006] 所述，ROI為感興趣區(qū)域；
[0007] 步驟二、對(duì)選定的初始片層組Sselected進(jìn)行優(yōu)化，得到優(yōu)化后片層集合Soptimized;
[000引步驟S、對(duì)優(yōu)化后片層集合Snptimized中的片層的ROI進(jìn)行闊值自動(dòng)選擇與半月板目 標(biāo)提??；
[0009]步驟四、根據(jù)提取出的半月板目標(biāo)巧在優(yōu)化后片層集合SDPtimized中進(jìn)行半月板 原始信息提取，得到包含分割出的半月板原始信息的ROI的片層Tk。
[001日]發(fā)明效果
[0011]實(shí)驗(yàn)中采用的=組實(shí)驗(yàn)數(shù)據(jù)是由哈爾濱醫(yī)科大學(xué)附屬醫(yī)院提供，采集于=位年齡 分布在30-50之間、于2014年在哈爾濱醫(yī)科大學(xué)附屬醫(yī)院進(jìn)行膝關(guān)節(jié)磁共振檢查的患者，圖 像為T(mén)l加權(quán)(TR= 500ms ,TE= 17ms)，片層間隔為4mm，像素尺寸為51巧512。為尊重用戶隱 私，成像序列分別稱為Sl、S2、S3，每個(gè)序列包含圖像36張。圖像經(jīng)過(guò)預(yù)處理已將DI COM格式 轉(zhuǎn)換為易于計(jì)算機(jī)處理的PNG格式。
[0012] 本實(shí)驗(yàn)的計(jì)算機(jī)配置為：cpu:Inter巧-4200M@2.5GHz，內(nèi)存：4G，操作系統(tǒng) WindowslO，處理軟件：Mat lab2013a O
[001引結(jié)果如圖I至圖5、表1、表2所示。圖I為利用灰度分布信息對(duì)圖像序列SI的均值圖 像進(jìn)行ROI定位的效果圖，添加了局部極小值捜索步驟，根據(jù)Matlab連通域提取函數(shù)獲得的 目標(biāo)0/k,m位置信息W及膝關(guān)節(jié)成像特點(diǎn)，可見(jiàn)此方法能準(zhǔn)確的定位半月板并且無(wú)丟失的提 取R0I。圖2為利用灰度分布信息對(duì)圖像序列的均值圖像進(jìn)行ROI定位的效果圖，未添加局部 極小值捜索步驟的結(jié)果對(duì)比圖，圖2a為未添加局部極小值捜索步驟時(shí)，利用灰度分布信息 對(duì)圖像序列Sl的均值圖像進(jìn)行ROI定位的效果圖，可見(jiàn)未添加局部極小值捜索步驟時(shí)ROI定 位會(huì)由于成像因素的影響(膝關(guān)節(jié)在圖像中的位置、膝關(guān)節(jié)彎曲程度等），出現(xiàn)在圖像兩端， 無(wú)法實(shí)現(xiàn)半月板準(zhǔn)確定位，圖化、圖2c分別為未添加局部極小值捜索步驟時(shí)，利用灰度分布 信息對(duì)圖像序列S2、S3的均值圖像進(jìn)行ROI定位的效果圖，可見(jiàn)運(yùn)種成像因素造成的影響并 不是總會(huì)出現(xiàn)，但出現(xiàn)幾率較高。圖3a為序列Sl中第5片層分割出的ROI圖像序列示例圖；圖 3b為序列Sl中第7片層分割出的ROI圖像序列示例圖；圖3c為序列Sl中第9片層分割出的ROI 圖像序列示例圖；圖3d為序列Sl中第11片層分割出的ROI圖像序列示例圖；圖3e為序列Sl中 第13片層分割出的ROI圖像序列示例圖；可見(jiàn)ROI的選擇可完整的包含半月板結(jié)構(gòu)，并且使 得面積盡可能減小便于運(yùn)算、降低誤提取概率。圖4a為序列Sl中第5片層半月板目標(biāo)提取結(jié) 果示例圖；圖4b為序列Sl中第7片層半月板目標(biāo)提取結(jié)果示例圖；圖4c為序列Sl中第9片層 半月板目標(biāo)提取結(jié)果示例圖；圖4d為序列Sl中第11片層半月板目標(biāo)提取結(jié)果示例圖；圖4e 為序列SI中第13片層半月板目標(biāo)提取結(jié)果示例圖；可見(jiàn)半月板結(jié)構(gòu)被完整提取出來(lái)，并且 已進(jìn)行闊值化處理和空桐填充，此時(shí)半月板目標(biāo)只有形狀信息(失去位置、灰度信息），從第 11片層結(jié)果可見(jiàn)存在誤提取現(xiàn)象（由于原目標(biāo)較小不易展示，部分圖像被放大）。圖5a為序 列Sl中第5片層半月板信息提取結(jié)果示例圖；圖化為序列Sl中第7片層半月板信息提取結(jié)果 示例圖；圖5c為序列Sl中第9片層半月板信息提取結(jié)果示例圖；圖5d為序列Sl中第11片層半 月板信息提取結(jié)果示例圖；圖5e為序列Sl中第13片層半月板信息提取結(jié)果示例圖；可見(jiàn)半 月板信息（包括形狀、位置關(guān)系、灰度紋理信息）被完整地從原MRI中提取出來(lái)，圖5曰、圖5b、 圖5c、圖5d、圖5e呈現(xiàn)的是整個(gè)ROI區(qū)域，從第11片層結(jié)果可見(jiàn)存在誤提取現(xiàn)象(該部分為灰 度相似，且滿足形狀約束的非半月板組織）。圖6a為添加了基于醫(yī)學(xué)先驗(yàn)知識(shí)結(jié)合位置關(guān)系 和目標(biāo)數(shù)量進(jìn)行濾除的步驟后，序列Sl中第11片層目標(biāo)誤提取現(xiàn)象消失對(duì)比圖；圖化為添 加了基于醫(yī)學(xué)先驗(yàn)知識(shí)結(jié)合位置關(guān)系和目標(biāo)數(shù)量進(jìn)行濾除的步驟后，序列Sl中第11片層目 標(biāo)誤提取現(xiàn)象消失結(jié)果圖。
[0014]表1半月板分割結(jié)果評(píng)估指標(biāo)
[0015]
[i
[i
[0018] ROI分割結(jié)果如圖I所示，與未添加局部極小值捜索的方法結(jié)果對(duì)比如圖2a-2c所 示，可見(jiàn)實(shí)現(xiàn)了存在兩端黑影干擾時(shí)的ROI準(zhǔn)確定位。半月板目標(biāo)提取結(jié)果如圖3a-3e所示， 與圖4a-4e(未添加誤提取目標(biāo)濾除）對(duì)比，可見(jiàn)增加的步驟可降低誤提取幾率。圖5a-f5e顯 示了最終半月板分割的結(jié)果。利用專家對(duì)MRI序列中半月板目標(biāo)的人工標(biāo)注作為參考標(biāo)準(zhǔn)， 計(jì)算出方法的性能指標(biāo):提取半月板目標(biāo)的準(zhǔn)確度、目標(biāo)提取的召回率，如表1所示。由于方 法中計(jì)算出的ROI面積W及闊值捜索范圍小于已有算法的1/2,此部分運(yùn)算效率為原有的二 倍。實(shí)驗(yàn)中還分析了對(duì)于不同形狀的半月板目標(biāo)的提取性能，結(jié)果如表2所示。
【附圖說(shuō)明】
[0019] 圖1為添加了局部極小值捜索步驟時(shí)，利用灰度分布信息對(duì)圖像序列Sl的均值圖 像進(jìn)行ROI定位的效果圖；
[0020] 圖2a為未添加局部極小值捜索步驟時(shí)，利用灰度分布信息對(duì)圖像序列Sl的均值圖 像進(jìn)行ROI定位的效果圖；
[0021] 圖化為未添加局部極小值捜索步驟時(shí)，利用灰度分布信息對(duì)圖像序列S2的均值圖 像進(jìn)行ROI定位的效果圖；
[0022] 圖2c為未添加局部極小值捜索步驟時(shí)，利用灰度分布信息對(duì)圖像序列S3的均值圖 像進(jìn)行ROI定位的效果圖；
[0023] 圖3a為序列Sl中第5片層分割出的ROI圖像序列示例圖；
[0024] 圖3b為序列Sl中第7片層分割出的ROI圖像序列示例圖；
[0025] 圖3c為序列Sl中第9片層分割出的ROI圖像序列示例圖；
[0026] 圖3d為序列Sl中第11片層分割出的ROI圖像序列示例圖；
[0027] 圖3e為序列Sl中第13片層分割出的ROI圖像序列示例圖；
[0028] 圖4a為序列Sl中第5片層半月板目標(biāo)提取結(jié)果示例圖；
[0029] 圖4b為序列Sl中第7片層半月板目標(biāo)提取結(jié)果示例圖；
[0030] 圖4c為序列Sl中第9片層半月板目標(biāo)提取結(jié)果示例圖；
[0031] 圖4d為序列Sl中第11片層半月板目標(biāo)提取結(jié)果示例圖；
[0032] 圖4e為序列Sl中第13片層半月板目標(biāo)提取結(jié)果示例圖；
[0033] 圖5a為序列Sl中第5片層半月板信息提取結(jié)果示例圖；
[0034] 圖化為序列Sl中第7片層半月板信息提取結(jié)果示例圖；
[0035] 圖5c為序列Sl中第9片層半月板信息提取結(jié)果示例圖；
[0036] 圖5d為序列Sl中第11片層半月板信息提取結(jié)果示例圖；
[0037] 圖5e為序列Sl中第13片層半月板信息提取結(jié)果示例圖；
[0038] 圖6a為圖4d經(jīng)過(guò)濾除誤提取目標(biāo)處理后的對(duì)比圖；
[0039] 圖化為圖5d經(jīng)過(guò)濾除誤提取目標(biāo)處理后的對(duì)比圖；
[0040] 圖7為子矩陣平均灰度曲線圖，Mean intensity for submatrix為子矩陣的灰度 平均值，Mean intensity為平均灰度，No.of submatrix為子矩陣的下標(biāo)；
[0041] 圖8為歸一化后的LSE的分布圖，橫坐標(biāo)是n縱坐標(biāo)是每個(gè)子矩陣計(jì)算的LSE;LSE between ROIs and mask為ROI和模板間的LSE，No.of ROIs為ROI的下標(biāo)；
[0042] 圖9為本發(fā)明流程圖；
[0043] 圖10為W均值化圖像Imean縱向中軸線為中屯、兩側(cè)各90像素寬度區(qū)域作為待處理 區(qū)域，沿待處理區(qū)域縱向中軸線垂直方向分割為16行X 180列像素大小的子矩陣示意圖。
【具體實(shí)施方式】
【具體實(shí)施方式】 [0044] 一:結(jié)合圖9說(shuō)明本實(shí)施方式，本實(shí)施方式的一種基于形狀約束的膝 關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取方法，具體是按照W下步驟制備的：
[0045] 可能包含半月板的片層可通過(guò)表3中信息初步選定，稱為Sselected,所有片層稱為 Sallo
[0046] 表3包含半月板的片層分布
[0047]
[0048] 半月板在膝關(guān)節(jié)MRI中出現(xiàn)的位置寬度，基本固定，通過(guò)經(jīng)驗(yàn)可確定其寬度可選擇 在MRI中，中線左右寬度180像素的區(qū)域。
[0049]根據(jù)經(jīng)驗(yàn)半月板高度一般在10-20像素之間。根據(jù)分析膝關(guān)節(jié)有可能出現(xiàn)的最大 傾斜，通過(guò)計(jì)算可完全包含半月板的最小可接受ROI為6府X18咧(計(jì)算方法:假設(shè)W-種極端 情況為例，兩半月板上下錯(cuò)開(kāi)50%厚度（10像素），平均圖像總高度達(dá)到30像素，再取極端情 況，定位小片定位ROI在一個(gè)半月板遠(yuǎn)離另一個(gè)的一半處，貝化OI-半寬度為30將完全包含 該半月板5像素的厚度，另一半30像素，也足夠包含另一邊1日+1日=30)。
[(K)加]步驟一、ROI (感興趣區(qū)域)定位，選定初始片層組Sselected;
[0051] 所述，ROI為感興趣區(qū)域；
[0052] 步驟二、對(duì)選定的初始片層組Sselected進(jìn)行優(yōu)化，得到優(yōu)化后片層集合Soptimized;
[0053] 步驟S、對(duì)優(yōu)化后片層集合SDPtimized中的片層的ROI進(jìn)行闊值自動(dòng)選擇與半月板目 標(biāo)提取；
[0054] 步驟四、根據(jù)提取出的半月板目標(biāo)巧r在優(yōu)化后片層集合Snptimized中進(jìn)行半月板 原始信息提取，得到包含分割出的半月板原始信息的ROI的片層化。
【具體實(shí)施方式】 [0055] 二:本實(shí)施方式與一不同的是:所述步驟一中ROI (感 興趣區(qū)域)定位，選定初始片層組Sselected ;具體過(guò)程為：
[0056] (1)計(jì)算選定的初始片層的均值化圖像，記為Imean;
[0057] 其中，所述選定為人為根據(jù)臨床解剖學(xué)知識(shí)選定;選定的初始片層為512*512像 素；
[00 日引 Imean - ITlGan { Ik( i , j ) I Ik E Sselected} (I)
[0059] 其中，Ik表示選定的初始片層組Sselected中的下標(biāo)為k的片層；i為1-512像素橫坐 標(biāo)，j為1-512像素縱坐標(biāo)；
[0060] (2) W均值化圖像Imean縱向中軸線為中屯、兩側(cè)各90像素寬度區(qū)域作為待處理區(qū) 域，沿待處理區(qū)域縱向中軸線垂直方向分割為16行X 180列像素大小的子矩陣，如圖10;
[0061] 第n個(gè)子矩陣表示為Pn(iJ)，用于定位半月板垂直位置，子矩陣高度的選擇接近 根據(jù)經(jīng)驗(yàn)的平均半月板厚度，同時(shí)與圖像整體尺寸有關(guān)(512巧12像素）；n取值為1-32;
[0062] (3)計(jì)算每個(gè)子矩陣的平均灰度值，如圖7,并繪制出橫坐標(biāo)為n，縱坐標(biāo)為灰度值 的子矩陣平均灰度曲線圖；
[0063] Np(n)=mean(Pn(i,j)),nG(l,32) (2)
[0064] 式中，Np(n)為第n個(gè)子矩陣的平均灰度值；
[0065] (4)根據(jù)經(jīng)驗(yàn)知識(shí)設(shè)置子矩陣的平均灰度值的局部極小值捜索范圍為[30%， 70%]，在范圍為[30% ,70% ]內(nèi)捜索平均灰度值最小值對(duì)應(yīng)的子矩陣作為包含半月板小 片；
[0066]
傅
[0067] 式中，1^為平均灰度值最小的子矩陣的下標(biāo)，11^￡11;
[006引將巧4/J)的縱向的高度/2+1行C巧.(z'J)的縱向的高度/2+1行為縱向的高度/2向 下移動(dòng)一行)作為ROI窗口的中屯、行，中屯、行向上選取30行，向下29行作為實(shí)際選定的ROI窗 口，將實(shí)際選定的ROI窗口作用于全部片層獲得ROI圖像序列礦''(;'.八取值為為1-180 像素橫坐標(biāo)，j取值為為1-60像素縱坐標(biāo)；
[0069]其中，/r (z'J)為下標(biāo)為k的片層的ROI ,Sail為全部片層集合。
[0070]其它步驟及參數(shù)與【具體實(shí)施方式】一相同。
[0071 ]【具體實(shí)施方式】=:本實(shí)施方式與【具體實(shí)施方式】一或二不同的是:所述步驟二中對(duì) 選定的初始片層組Sseleeted進(jìn)行優(yōu)化，得到優(yōu)化后片層集合Snptimized ;具體過(guò)程為：
[0072] (1)計(jì)算選定的初始片層組SseleGted的ROI的均值作為模板& (i，j);
[007；3]
(4)
[0074] 其中，M為選定的初始片層組Sselected中的片層數(shù)；(Hi, j)為選定的初始片層組 Sselected的ROI的均值；（2)利用模板(I)QJ)匹配所有選定片層Ik的ROI，計(jì)算最小均方誤差 LSE，并進(jìn)行歸一化處理，如圖8;
[0075] 樹(shù)
[0076] ㈱
[0077] 將歸一化處理后LSE^化)〉60%的片層視為不包含半月板片層，剔除運(yùn)些片層，余 下的片層作為優(yōu)化后的片層集合SDPtimized。
[0078] 如果要選擇出典型的包含類S角形半月板的MRI，則可將條件變嚴(yán)格，設(shè)置LSE^ 化)〉40%或者更小，則靠近兩側(cè)和中部的片層將被刪除。
[0079] 其它步驟及參數(shù)與【具體實(shí)施方式】一或二相同。
【具體實(shí)施方式】 [0080] 四：本實(shí)施方式與一至=之一不同的是：步驟=中對(duì) 優(yōu)化后片層集合SDPtimized中的片層的ROI進(jìn)行闊值自動(dòng)選擇與半月板目標(biāo)提取;具體過(guò)程 為：
[0081] (1)計(jì)算優(yōu)化后片層集合Soptimized中的片層的ROI的均值圖像的平均灰度，作為闊 值選擇的參考值；
[0082]
(7)
[0083] 式中，n為優(yōu)化后片層集合Soptimized中的片層的ROI的均值圖像的平均灰度；
[0084] (2)選取闊值T = O.0511，對(duì)優(yōu)化后片層集合Soptimized中的片層的ROI進(jìn)行二值化處 理；
[0085]
(S)
[0086] 式中，化(i，j)為二值化處理后的圖像；
[0087] (3)將二值化處理后的圖像Bk(i J )的ROI中所有滿足八連通的獨(dú)立目標(biāo)利用 Matlab連通域提取函數(shù)全部提取出來(lái)，進(jìn)行濾波約束處理，記為目標(biāo)Ok,m，k為片層編號(hào)，m為 每片層中目標(biāo)編號(hào);此時(shí)目標(biāo)包含位置坐標(biāo)；
[0088] (4)首先進(jìn)行過(guò)大或過(guò)小目標(biāo)的濾除，將面積不滿足表4中面積約束條件的目標(biāo)濾 除，初步減少噪聲像素的干擾。
[0089] 當(dāng)目標(biāo)0k,m不滿足面積為Area(0k,m)〉50像素，Area(0k,m)<800像素時(shí)，濾除目標(biāo) Ok, m;
[0090] 表4目標(biāo)形狀約束條件
[0091]
[0092] (5)針對(duì)濾除目標(biāo)Ok,m后每一個(gè)剩余目標(biāo)，計(jì)算剩余目標(biāo)的軸比和邊界盒比，軸比 為主軸/副軸，邊界盒比為邊界長(zhǎng)度/寬度，濾除不滿足軸比和邊界盒比的目標(biāo)，得到滿足半 月板形狀約束條件的目標(biāo)O^k,"，每個(gè)片層上ROI中目標(biāo)數(shù)量記為姑，,，下標(biāo)t為闊值變化標(biāo) 記，1《t《19;形狀約束條件見(jiàn)表4。
[0093] 其中，軸比的目標(biāo)為主軸/副軸〉1，主軸/副軸<6,邊界盒比目標(biāo)為長(zhǎng)度/寬度〉1，長(zhǎng) 度/寬度巧;半月板形狀約束條件為軸比目標(biāo)為主軸/副軸〉1，主軸/副軸<6,邊界盒比目標(biāo) 為長(zhǎng)度/寬度〉1，長(zhǎng)度/寬度巧；
[0094]
(9)
[00M] (6)根據(jù)Matlab連通域提取函數(shù)獲得的目標(biāo)OYm位置信息W及膝關(guān)節(jié)成像特點(diǎn)， 濾除誤提取的目標(biāo)，得到濾除誤提取的目標(biāo)后的提取結(jié)果〇"k,m(若出現(xiàn)剩余目標(biāo)超過(guò)3個(gè)， 則判定兩側(cè)目標(biāo)為半月板，其余目標(biāo)為關(guān)節(jié)軟骨殘余部分）；
[0096]
(10)
[0097] 式中，0"k,"為濾除誤提取的目標(biāo)后的提取結(jié)果；
[009引（7)記錄每個(gè)片層上ROI中提取的半月板目標(biāo)數(shù)量qk, t;
[0099] qk,t = no.of (0"k'm) (11)
[0100] W及闊值T對(duì)應(yīng)的ROI序列總共提取出的半月板目標(biāo)數(shù)量；
[0101]
(12)
[0102] (8)改變闊值，重復(fù)步驟(3)-(7)，記錄用每個(gè)闊值處理時(shí)提取出的半月板目標(biāo)總 數(shù)量qt，闊值范圍為T(mén) = 0.0時(shí)~0.4〇ri，每次增長(zhǎng)0.OSn，其中下標(biāo)t = 1，2，......,8；
[0103] (9)選擇使提取出的半月板目標(biāo)總數(shù)量qt最多的闊值作為自適應(yīng)闊值選擇的最優(yōu) 結(jié)果；
[0104]
(13)
[0105] 乍為二值化處理的最優(yōu)闊值，T氣4應(yīng)的提取結(jié)果〇"k,m即為最優(yōu)半月板目標(biāo)提 取結(jié)果，此時(shí)提取半月板目標(biāo)記為巧:f。
[0106] 其它步驟及參數(shù)與【具體實(shí)施方式】一至=之一相同。
【具體實(shí)施方式】 [0107] 五:本實(shí)施方式與一至四之一不同的是:所述T巧X值 為為0.1 On~0.3011。通過(guò)實(shí)驗(yàn)我們發(fā)現(xiàn)對(duì)于我們的Tl加權(quán)MRI數(shù)據(jù)，最優(yōu)T巧別直為0. l〇ri~ 0.3〇ri。
[0108] 其它步驟及參數(shù)與【具體實(shí)施方式】一至四之一相同。
【具體實(shí)施方式】 [0109] 六:本實(shí)施方式與一至五之一不同的是:所述步驟四 中根據(jù)提取出的半月板目標(biāo)0品'在優(yōu)化后片層集合SDPtimized中進(jìn)行半月板原始信息提取，得 到包含分割出的半月板原始信息的ROI的片層化;具體過(guò)程為：
[0110] (1)根據(jù)半月板目標(biāo)《r，計(jì)算0汾的凸包巧:f?'(填補(bǔ)闊值處理可能產(chǎn)生的目標(biāo)內(nèi) 部空桐)作為模板。
[0111] (2)利用第二階段步驟中Matlab連通域提取函數(shù)記錄的每一個(gè)化r的位置坐標(biāo)，獲 取半月板目標(biāo)對(duì)應(yīng)到原片層中的位置坐標(biāo)。
[0112] (3)制作與ROI區(qū)域大小(60行*180列像素)相同的模板，其中對(duì)應(yīng)于凸化處理后的 半月板目標(biāo)化T*"'"的區(qū)域?yàn)榭眨渌麉^(qū)域灰度設(shè)置為255(白色，便于展示與觀察）；
[0113]
(14)
[0114] 化即為最終僅包含分割出的半月板原始信息的ROI的片層。
[0115] 其它步驟及參數(shù)與【具體實(shí)施方式】一至五之一相同。
[0116] 采用W下實(shí)施例驗(yàn)證本發(fā)明的有益效果：
[0117] 實(shí)施例一；
[0118] 本實(shí)施例一種基于形狀約束的膝關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取方法具體 是按照W下步驟制備的：
[0119] 實(shí)驗(yàn)中采用的=組實(shí)驗(yàn)數(shù)據(jù)是由哈爾濱醫(yī)科大學(xué)附屬醫(yī)院提供，采集于=位年齡 分布在30-50之間、于2014年在哈爾濱醫(yī)科大學(xué)附屬醫(yī)院進(jìn)行膝關(guān)節(jié)磁共振檢查的患者，圖 像為T(mén)l加權(quán)(TR=SOOms ,TE= 17ms)，片層間隔為4mm，像素尺寸為51巧512。為尊重用戶隱 私，成像序列分別稱為51、52、53，每個(gè)序列包含圖像36張。圖像經(jīng)過(guò)預(yù)處理已將01〇)1格式 轉(zhuǎn)換為易于計(jì)算機(jī)處理的PNG格式。
[0120] 本實(shí)驗(yàn)的計(jì)算機(jī)配置為：cpu:Inter巧-4200M@2.5GHz，內(nèi)存：4G，操作系統(tǒng) WindowslO，處理軟件：Matlab2013a。
[0121] 結(jié)果如圖1至圖5、表1、表2所示。圖1為利用灰度分布信息對(duì)圖像序列SI的均值圖 像進(jìn)行ROI定位的效果圖，添加了局部極小值捜索步驟，根據(jù)Matlab連通域提取函數(shù)獲得的 目標(biāo)〇/k,m位置信息W及膝關(guān)節(jié)成像特點(diǎn)，可見(jiàn)此方法能準(zhǔn)確的定位半月板并且無(wú)丟失的提 取R0I。圖2為利用灰度分布信息對(duì)圖像序列的均值圖像進(jìn)行ROI定位的效果圖，未添加局部 極小值捜索步驟的結(jié)果對(duì)比圖，圖2a為未添加局部極小值捜索步驟時(shí)，利用灰度分布信息 對(duì)圖像序列Sl的均值圖像進(jìn)行ROI定位的效果圖，可見(jiàn)未添加局部極小值捜索步驟時(shí)ROI定 位會(huì)由于成像因素的影響(膝關(guān)節(jié)在圖像中的位置、膝關(guān)節(jié)彎曲程度等），出現(xiàn)在圖像兩端， 無(wú)法實(shí)現(xiàn)半月板準(zhǔn)確定位，圖化、圖2c分別為未添加局部極小值捜索步驟時(shí)，利用灰度分布 信息對(duì)圖像序列S2、S3的均值圖像進(jìn)行ROI定位的效果圖，可見(jiàn)運(yùn)種成像因素造成的影響并 不是總會(huì)出現(xiàn)，但出現(xiàn)幾率較高。圖3a為序列Sl中第5片層分割出的ROI圖像序列示例圖；圖 3b為序列Sl中第7片層分割出的ROI圖像序列示例圖；圖3c為序列Sl中第9片層分割出的ROI 圖像序列示例圖；圖3d為序列Sl中第11片層分割出的ROI圖像序列示例圖；圖3e為序列Sl中 第13片層分割出的ROI圖像序列示例圖；可見(jiàn)ROI的選擇可完整的包含半月板結(jié)構(gòu)，并且使 得面積盡可能減小便于運(yùn)算、降低誤提取概率。圖4a為序列Sl中第5片層半月板目標(biāo)提取結(jié) 果示例圖；圖4b為序列SI中第7片層半月板目標(biāo)提取結(jié)果示例圖；圖4c為序列SI中第9片層 半月板目標(biāo)提取結(jié)果示例圖；圖4d為序列Sl中第11片層半月板目標(biāo)提取結(jié)果示例圖；圖4e 為序列SI中第13片層半月板目標(biāo)提取結(jié)果示例圖；可見(jiàn)半月板結(jié)構(gòu)被完整提取出來(lái)，并且 已進(jìn)行闊值化處理和空桐填充，此時(shí)半月板目標(biāo)只有形狀信息(失去位置、灰度信息），從第 11片層結(jié)果可見(jiàn)存在誤提取現(xiàn)象（由于原目標(biāo)較小不易展示，部分圖像被放大）。圖5a為序 列Sl中第5片層半月板信息提取結(jié)果示例圖；圖化為序列Sl中第7片層半月板信息提取結(jié)果 示例圖；圖5c為序列Sl中第9片層半月板信息提取結(jié)果示例圖；圖5d為序列Sl中第11片層半 月板信息提取結(jié)果示例圖；圖5e為序列Sl中第13片層半月板信息提取結(jié)果示例圖；可見(jiàn)半 月板信息（包括形狀、位置關(guān)系、灰度紋理信息）被完整地從原MRI中提取出來(lái)，圖5曰、圖5b、 圖5c、圖5d、圖5e呈現(xiàn)的是整個(gè)ROI區(qū)域，從第11片層結(jié)果可見(jiàn)存在誤提取現(xiàn)象(該部分為灰 度相似，且滿足形狀約束的非半月板組織）。圖6a為添加了基于醫(yī)學(xué)先驗(yàn)知識(shí)結(jié)合位置關(guān)系 和目標(biāo)數(shù)量進(jìn)行濾除的步驟后，序列Sl中第11片層目標(biāo)誤提取現(xiàn)象消失對(duì)比圖；圖化為添 加了基于醫(yī)學(xué)先驗(yàn)知識(shí)結(jié)合位置關(guān)系和目標(biāo)數(shù)量進(jìn)行濾除的步驟后，序列Sl中第11片層目 標(biāo)誤提取現(xiàn)象消失結(jié)果圖。
[0122] 表1半月板分割結(jié)果評(píng)估指標(biāo)
[0123]
[0124] 表2不同形態(tài)半月板檢測(cè)性能比較
[0125]
[0126] ROI分割結(jié)果如圖I所示，與未添加局部極小值捜索的方法結(jié)果對(duì)比如圖2所示，可 見(jiàn)實(shí)現(xiàn)了存在兩端黑影干擾時(shí)的ROI準(zhǔn)確定位。半月板目標(biāo)提取結(jié)果如圖3a-3e所示，與圖 4a-4e(未添加誤提取目標(biāo)濾除）對(duì)比，可見(jiàn)增加的步驟可降低誤提取幾率。圖5曰-扣顯示了 最終半月板分割的結(jié)果。利用專家對(duì)MRI序列中半月板目標(biāo)的人工標(biāo)注作為參考標(biāo)準(zhǔn)，計(jì)算 出方法的性能指標(biāo):提取半月板目標(biāo)的準(zhǔn)確度、目標(biāo)提取的召回率，如表1所示。由于方法中 計(jì)算出的ROI面積W及闊值捜索范圍小于已有算法的1/2,此部分運(yùn)算效率為原有的二倍。 實(shí)驗(yàn)中還分析了對(duì)于不同形狀的半月板目標(biāo)的提取性能，結(jié)果如表2所示。
[0127] 本發(fā)明還可有其它多種實(shí)施例，在不背離本發(fā)明精神及其實(shí)質(zhì)的情況下，本領(lǐng)域 技術(shù)人員當(dāng)可根據(jù)本發(fā)明作出各種相應(yīng)的改變和變形，但運(yùn)些相應(yīng)的改變和變形都應(yīng)屬于 本發(fā)明所附的權(quán)利要求的保護(hù)范圍。
【主權(quán)項(xiàng)】
1. 一種基于形狀約束的膝關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取方法，其特征在于具體 是按照以下步驟進(jìn)行的： 步驟一、ROI定位，選定初始片層組Sselected ; 所述，ROI為感興趣區(qū)域； 步驟二、對(duì)選定的初始片層組&^。^進(jìn)行優(yōu)化，得到優(yōu)化后片層集合SoptWmd; 步驟三、對(duì)優(yōu)化后片層集合SoptWzed中的片層的ROI進(jìn)行閾值自動(dòng)選擇與半月板目標(biāo)提 ??； 步驟四、根據(jù)提取出的半月板目標(biāo)0ΗΓ在優(yōu)化后片層集合中進(jìn)行半月板原始 信息提取，得到包含分割出的半月板原始信息的R0I的片層Tk。2. 根據(jù)權(quán)利要求1所述一種基于形狀約束的膝關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取方 法，其特征在于:所述步驟一中R0I定位，選定初始片層組S selec；ted;具體過(guò)程為： (1) 計(jì)算選定的初始片層的均值化圖像，記為Im_; 選定的初始片層為512*512像素； Imean - ??Θ??Ι {lk(i,j)|lkG Sselected } ( 1 ) 其中，Ik表示選定的初始片層組Sselected中的下標(biāo)為k的片層；i為1-512像素橫坐標(biāo)，j為 1-512像素縱坐標(biāo)； (2) 以均值化圖像Imean縱向中軸線為中心兩側(cè)各90像素寬度區(qū)域作為待處理區(qū)域，沿待 處理區(qū)域縱向中軸線垂直方向分割為16行X 180列像素大小的子矩陣； 第η個(gè)子矩陣表示SPn( i，j);η取值為1-32; (3) 計(jì)算每個(gè)子矩陣的平均灰度值， Np(n) =mean(Pn(i , j)) ,ne (1,32) (2) 式中，NP(n)為第n個(gè)子矩陣的平均灰度值； (4) 設(shè)置子矩陣的平均灰度值的局部極小值搜索范圍為[30%，70%]，在范圍為[30%， 70% ]內(nèi)搜索平均灰度值最小值對(duì)應(yīng)的子矩陣作為包含半月板小片；(3) 式中，為平均灰度值最小的子矩陣的下標(biāo)， 將gU的縱向的高度/2+1行作為ROI窗口的中心行，中心行向上選取30行，向下29行 作為實(shí)際選定的R0I窗口，將實(shí)際選定的R0I窗口作用于全部片層獲得R0I圖像序列 Saii，i取值為為1_180像素橫坐標(biāo)，j取值為為1_60像素縱坐標(biāo)； 其中，(/,/)為下標(biāo)為k的片層的ROI，Sall為全部片層集合。3. 根據(jù)權(quán)利要求2所述一種基于形狀約束的膝關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取方 法，其特征在于:所述步驟二中對(duì)選定的初始片層組&^。^進(jìn)行優(yōu)化，得到優(yōu)化后片層集合 Soptimized ;具體過(guò)程為： (1 )計(jì)算選定的初始片層組Sselected的R0I的均值作為模板Φ ( i，j ); 其中，Μ為選定的初始片層組Sselecrted中的片層數(shù)； Φ (i，j)為選定的初始片層組Sselected的ROI的均值； (2)利用模板匹配所有選定片層Ik的ROI，計(jì)算最小均方誤差LSE，并進(jìn)行歸一化 處理； k ' " · ·將歸一化處理后LSElkpeo%的片層視為不包含半月板片層，剔除這些片層，余下的 片層作為優(yōu)化后的片層集合Soptimized。4.根據(jù)權(quán)利要求3所述一種基于形狀約束的膝關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取方 法，其特征在于:步驟三中對(duì)優(yōu)化后片層集合SoptWzed中的片層的R0I進(jìn)行閾值自動(dòng)選擇與 半月板目標(biāo)提取;具體過(guò)程為： (1) 計(jì)算優(yōu)化后片層集合中的片層的R0I的均值圖像的平均灰度，作為閾值選 擇的參考值；式中， η為優(yōu)化后片層集合Soptimized中的片層的R01的均值圖像的平均灰度； (2) 選取閾值τ = 0.05η，對(duì)優(yōu)化后片層集合Soptimized中的片層的R0I進(jìn)行二值化處理；式中， Bk(i，j)為二值化處理后的圖像； (3) 將二值化處理后的圖像Bk(i，j)的R0I中所有滿足八連通的獨(dú)立目標(biāo)利用Matlab連 通域提取函數(shù)全部提取出來(lái)，進(jìn)行濾波約束處理，記為目標(biāo)Ok, m，k為片層編號(hào)，m為每片層中 目標(biāo)編號(hào)； (4) 當(dāng)目標(biāo)Ok,m不滿足面積為Area(0k,m)>50像素，Area(0k,m)〈800像素時(shí)，濾除目標(biāo)Ok,m; (5) 針對(duì)濾除目標(biāo)0k,m后每一個(gè)剩余目標(biāo)，計(jì)算剩余目標(biāo)的軸比和邊界盒比，軸比為主 軸/副軸，邊界盒比為邊界長(zhǎng)度/寬度，濾除不滿足軸比和邊界盒比的目標(biāo)，得到滿足半月板 形狀約束條件的目標(biāo)(/k, m，每個(gè)片層上R0I中目標(biāo)數(shù)量記為ft，下標(biāo)t為閾值變化標(biāo)記，K t^l9； 其中，軸比的目標(biāo)為主軸/副軸>1，主軸/副軸〈6,邊界盒比目標(biāo)為長(zhǎng)度/寬度>1，長(zhǎng)度/ 寬度〈5;半月板形狀約束條件為軸比目標(biāo)為主軸/副軸>1，主軸/副軸〈6,邊界盒比目標(biāo)為長(zhǎng) 度/寬度>1，長(zhǎng)度/寬度〈5;(6) 根據(jù)Matlab連通域提取函數(shù)獲得的目標(biāo)亇k,m位置信息以及膝關(guān)節(jié)成像特點(diǎn)，濾除 誤提取的目標(biāo)，得到濾除誤提取的目標(biāo)后的提取結(jié)果〇〃k, m; (10) 式中，〇\,?為濾除誤提取的目標(biāo)后的提取結(jié)果； (7) 記錄每個(gè)片層上ROI中提取的半月板目標(biāo)數(shù)量qk,t; qk，t = no.of(0〃k，m) (11) 以及閾值τ對(duì)應(yīng)的ROI序列總共提取出的半月板目標(biāo)數(shù)量；(8) 改變閾值，重復(fù)步驟(3) - (7)，記錄用每個(gè)閾值處理時(shí)提取出的半月板目標(biāo)總數(shù)量 qt，閾值范圍為τ = 0.05η~〇.40η，每次增長(zhǎng)〇.〇5n，其中下標(biāo)t = 1，2，......,8； (9) 選擇使提取出的半月板目標(biāo)總數(shù)量qt最多的閾值作為自適應(yīng)閾值選擇的最優(yōu)結(jié)果；以f作為二值化處理的最優(yōu)閾值，^對(duì)應(yīng)的提取結(jié)果〇〃k,m即為最優(yōu)半月板目標(biāo)提取結(jié) 果，此時(shí)提取半月板目標(biāo)記為。5. 根據(jù)權(quán)利要求4所述一種基于形狀約束的膝關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取方 法，其特征在于:所述f取值為為ο. 10η~〇. 30η。6. 根據(jù)權(quán)利要求5所述一種基于形狀約束的膝關(guān)節(jié)磁共振圖像序列半月板自動(dòng)提取方 法，其特征在于:所述步驟四中根據(jù)提取出的半月板目標(biāo)"Γ/在優(yōu)化后片層集合中 進(jìn)行半月板原始信息提取，得到包含分割出的半月板原始信息的R0I的片層Tk;具體過(guò)程 為： (1) 根據(jù)半月板目標(biāo)，計(jì)算0=的凸包〇；7'_1乍為模板； (2) 利用Matlab連通域提取函數(shù)記錄的每一個(gè)的位置坐標(biāo)，獲取半月板目標(biāo)對(duì)應(yīng)到 原片層中的位置坐標(biāo)； (3) 制作與roi區(qū)域大小相同的模板，其中對(duì)應(yīng)于凸化處理后的半月板目標(biāo)〇r!r?nvex的 區(qū)域?yàn)榭?，其他區(qū)域灰度設(shè)置為255;Tk即為最終僅包含分割出的半月板原始信息的R0I的片層。
【文檔編號(hào)】G06K9/34GK105956587SQ201610247856
【公開(kāi)日】2016年9月21日
【申請(qǐng)日】2016年4月20日
【發(fā)明人】李楊, 季云飛, 張寧, 齊月賓
【申請(qǐng)人】哈爾濱工業(yè)大學(xué)