面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù)處理方法
【專利摘要】本發(fā)明公開了一種面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù)處理方法，在圖形處理器中進(jìn)行以下步驟：步驟一、采集結(jié)腸部位的三維圖像數(shù)據(jù)，并對(duì)所述圖像數(shù)據(jù)進(jìn)行圖像預(yù)處理，得到結(jié)腸壁的擬合等值面曲面，在所述擬合等值面上選取體現(xiàn)結(jié)腸息肉特征的W個(gè)種子點(diǎn)；步驟二、逐步移動(dòng)所述種子點(diǎn)并將移動(dòng)后的種子點(diǎn)影到圖像真實(shí)等值面上，并形成W條曲率線；步驟三、將每一條所述曲率線進(jìn)行散播，生成曲率線集合；步驟四、在所述曲率線集合中篩選出體現(xiàn)結(jié)腸息肉形狀特征的特征曲率線，并將所述特征曲率線繪制到所述擬合等值面上，以突顯結(jié)腸息肉。本發(fā)明解決了現(xiàn)有技術(shù)中無法精確描述息肉的整體三維形狀特征的技術(shù)問題。
【專利說明】
面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù)處理方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明涉及無創(chuàng)檢測(cè)技術(shù)領(lǐng)域，特別是一種面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù) 處理方法。
【背景技術(shù)】
[0002] 虛擬結(jié)腸鏡是一種無創(chuàng)檢測(cè)技術(shù)。但是，基于CT或MR圖像的虛擬結(jié)腸鏡檢查十分 耗時(shí)，這就要求使用計(jì)算機(jī)輔助檢測(cè)系統(tǒng)對(duì)結(jié)腸圖像數(shù)據(jù)進(jìn)行預(yù)處理，以減少檢查時(shí)間?，F(xiàn) 有的結(jié)腸數(shù)據(jù)預(yù)處理方法主要基于圖像的標(biāo)量曲率，只能反映極小區(qū)域內(nèi)的三維形狀，無 法精確描述息肉的整體三維狀特征。并且，標(biāo)量曲率對(duì)數(shù)據(jù)中的噪聲比較敏感，導(dǎo)致最終的 檢測(cè)結(jié)果假陽性率很高。

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

[0003] 針對(duì)上述技術(shù)問題，本發(fā)明中提出了一種面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù)處 理方法，該方法使用基于曲率線的方法處理結(jié)腸息肉數(shù)據(jù)，能夠?qū)^大尺寸范圍內(nèi)的三維 形狀進(jìn)行分析，對(duì)結(jié)腸息肉的形狀特征的概括能力更強(qiáng)。同時(shí)，本發(fā)明描述的方法，還具有 穩(wěn)定性高、假陽性率低、自動(dòng)化程度高和適于并行化處理等優(yōu)點(diǎn)。本發(fā)明解決了現(xiàn)有技術(shù)中 無法精確描述息肉的整體三維狀特征的技術(shù)問題。
[0004] 為了實(shí)現(xiàn)根據(jù)本發(fā)明的這些目的和其它優(yōu)點(diǎn)，提供了一種面向虛擬結(jié)腸鏡的結(jié)腸 息肉圖像數(shù)據(jù)處理方法，包括：
[0005] 在圖形處理器中進(jìn)行以下步驟：
[0006] 步驟一、采集結(jié)腸部位的三維圖像數(shù)據(jù)，并對(duì)所述圖像數(shù)據(jù)進(jìn)行圖像預(yù)處理，得到 結(jié)腸壁的擬合等值面曲面，在所述擬合等值面上選取體現(xiàn)結(jié)腸息肉特征的W個(gè)種子點(diǎn)；
[0007] 步驟二、逐步移動(dòng)所述種子點(diǎn)并將移動(dòng)后的種子點(diǎn)影到圖像真實(shí)等值面上，連接 種子點(diǎn)及其在圖像真實(shí)等值面上的投影點(diǎn)形成W條曲率線；
[0008] 步驟三、將每一條所述曲率線進(jìn)行散播，生成曲率線集合；
[0009] 步驟四、在所述曲率線集合中篩選出體現(xiàn)結(jié)腸息肉形狀特征的特征曲率線，并將 所述特征曲率線繪制到所述擬合等值面上，以突顯結(jié)腸息肉。
[0010] 優(yōu)選的，所述步驟一中，圖形處理器讀入結(jié)腸部位的CT三維圖像f，分割出結(jié)腸壁 圖像，使用Marching Cubes算法計(jì)算生成分布結(jié)腸壁的所述擬合等值面曲面。
[0011]優(yōu)選的，所述步驟一中，將所述擬合等值面曲面映射到所述三維圖像的坐標(biāo)系中， 所述擬合等值面曲面由若干個(gè)頂點(diǎn)和邊組成，將頂點(diǎn)坐標(biāo)和頂點(diǎn)序號(hào)一一對(duì)應(yīng)。
[0012] 優(yōu)選的，所述步驟一中還包括:計(jì)算每個(gè)頂點(diǎn)的最大主曲率mi和最小主曲率m2，當(dāng) 每個(gè)頂點(diǎn)的豐曲率滿足以下備件時(shí)，該頂點(diǎn)作為候選種子點(diǎn)s ' ：
[0013]
[0014]其中，mm_是該頂點(diǎn)的平均曲率，mm_ = (mi+m2) /2，mt是預(yù)先選取的閾值，計(jì)算候選 種子點(diǎn)s'在其鄰域半徑r范圍內(nèi)所有相鄰頂點(diǎn)sx的主曲率，若sx的最大主曲率和最小主曲率 符號(hào)相反，則將候選種子點(diǎn)s'刪除，若Sx的最大主曲率和最小主曲率符號(hào)相同，則將該候選 種子點(diǎn)s'作為種子點(diǎn)s放入堆棧S中。
[0015] 優(yōu)選的，所述步驟二中包括以下步驟：
[0016] (I)、計(jì)算每個(gè)種子點(diǎn)s的最大主曲率ki和最小主曲率k2;
[0017] (II)、在每一個(gè)種子點(diǎn)s對(duì)應(yīng)的真實(shí)等值面的切平面上，沿著最大曲率方向逐步移 動(dòng)，其移動(dòng)步長(zhǎng)為L(zhǎng)，其定義為
其中，lu(i = l，2)是真實(shí)等值面的主曲 率，eL是設(shè)定的閾值；
[0018] (III)、對(duì)于S中的種子點(diǎn)s，其初始位置為P，該位置處對(duì)應(yīng)的體素灰度值為Ul，由 圖像處理器線程來計(jì)算該位置處的梯度向量I，若Ul小于真實(shí)等值面的體素灰度值U2,則定 義投影向量D為I，若Ul大于U2,則:
[0019] (VI)、將s在切平面上沿著最大主曲率移動(dòng)后的位置記作P⑴，P⑴沿著D的方向逐 步移動(dòng)，每一步移動(dòng)的步長(zhǎng)為d，其中，d定義為d = Lmax/20，在每一步移動(dòng)之后，重新獲取當(dāng) 前位置Sm的體素灰度，并記錄移動(dòng)前后的體素灰度分別是i前和i后；
[0020] (V)、若（itrU2)X(i后-U2K0,則停止移動(dòng)s，并計(jì)算sm在真實(shí)等值面上的投影位置 P，，
[0021]
[0022] 將P'放入堆棧S中，若種子點(diǎn)在切平面上移動(dòng)的補(bǔ)償LSLmax/lOO,或使用了四階 Runge-Kutta方法追蹤曲率線，則種子點(diǎn)在真實(shí)等值面上的投影過程可以忽略，在曲率線追 蹤停止之前，將投影位置P '作為種子點(diǎn)不斷放入堆棧S中，且將投影位置P '作為種子點(diǎn)回到 步驟(1)，連接在每個(gè)種子點(diǎn)s及其相應(yīng)的若干個(gè)投影位置P'形成一條曲率線；
[0023] (IV)、當(dāng)移動(dòng)后的種子點(diǎn)和其余的曲率線距離小于預(yù)設(shè)的閾值，或者種子點(diǎn)移動(dòng) 次數(shù)達(dá)到了預(yù)設(shè)的上限，則曲率線終止，將堆棧S中的種子點(diǎn)信息存儲(chǔ)到圖像處理器的共享 存儲(chǔ)器或顯存中。
[0024]優(yōu)選的，所述步驟三中，曲率線散播的具體步驟包括：
[0025]步驟A、在已經(jīng)生成的曲率線集合中任意選取一條曲率線1，選取1上的一個(gè)種子點(diǎn) s，真實(shí)等值面曲面的切線方向記作T，由g和T定義一個(gè)法平面N;
[0026] 步驟B、定義兩個(gè)向量VjPV2 ,V1 = TXg ,V1 = gXT，并開辟兩個(gè)圖像處理器線程"和 t2，tl和t2分別在Vl和V2的方向上分別放置兩個(gè)新增種子點(diǎn)Sl和S2，S1和S2到S的距離為ds，
[0027]
[0028] 若81到已有種子點(diǎn)的距離大于預(yù)設(shè)閾值τ，則將81放入堆棧S;若82到已有種子點(diǎn)的 距離大于預(yù)設(shè)閾值τ，則將 82放入堆棧S;不滿足條件的新增種子點(diǎn)對(duì)應(yīng)的線程終止計(jì)算，其 中，ki(i = l，2)是主曲率，ει是設(shè)定值；
[0029] 步驟C、分別以81和82為起始點(diǎn)，重復(fù)步驟一和二，生成新的曲率線UP12，UP1:^ 終止條件為:移動(dòng)后的種子點(diǎn)和其余的曲率線距離小于預(yù)設(shè)的閾值vt，或者種子點(diǎn)移動(dòng)次 數(shù)達(dá)到了預(yù)設(shè)的上限，其中，Vt為
[0030]
[0031]參量ε2是關(guān)于ε!的函數(shù)，可表達(dá)為 [0032]
[0033] 步驟D、曲率線UPl2生成后，刪除線程t#Pt2,對(duì)于UPl2,在UPl 2上分別選取一 個(gè)種子點(diǎn)，并重復(fù)步驟B和C，生成對(duì)應(yīng)曲率線;使用二叉樹結(jié)構(gòu)存儲(chǔ)種子點(diǎn)，同一深度的二 叉樹節(jié)點(diǎn)中的種子點(diǎn)在圖像處理器中是并行計(jì)算的；
[0034] 步驟E、重復(fù)上述步驟D，直至所生成的全部候選種子點(diǎn)到已有種子點(diǎn)的距離小于 預(yù)設(shè)閾值τ，曲率線散播結(jié)束，以上步驟全部完成后，生成曲率線集合Ln;然后使用最大曲 率，重復(fù)以上全部步驟，生成曲率線集合L m，Lm是所有沿著最大曲率方向生成的曲率線集 合，Ln是所有沿著最小曲率方向生成的曲率線集合，匕和1^為所有曲率線的集合。
[0035] 優(yōu)選的，所述步驟四中，篩選特征曲率線的步驟包括：
[0036] 步驟I、構(gòu)造雙曲面百分比參數(shù)HP，W5 = 100%，其中，nh是位于雙曲面上的曲 率線1的種子點(diǎn)，N為曲率線1上的全部種子點(diǎn)，若曲率線1上的種子點(diǎn)位于某一個(gè)潛在的雙 曲面，則其主曲率的符號(hào)相反；
[0037]步驟J、使用Winding angle表征曲率線的閉合性，記種子點(diǎn)s在移動(dòng)前后分別記作 slPs+，經(jīng)過s的切平面記作T，向量s-s在T上的投影向量記作（si) '，向量ss+在T上的投影 向量記作（Ss+) '，Winding angle是投影向量(si) '到（SS+) '的夾角，對(duì)篩選后的每條曲率 線計(jì)算其Winding angle，選擇Winding angle最大的若干條曲率線，并且進(jìn)一步要求曲率 線的起點(diǎn)和終點(diǎn)小于預(yù)設(shè)的距離閾值；
[0038]步驟K、計(jì)算篩選后的曲率線的平均半徑，其定義為曲率線平均中心到曲率線上各 種子點(diǎn)的平均距離，取平均曲率大于預(yù)設(shè)平均曲率閾值rt的曲率線作為特征曲率線。
[0039]本發(fā)明至少包括以下有益效果：
[0040] 1、本發(fā)明提出的方法，使用曲率線這一向量形式對(duì)圖像數(shù)據(jù)進(jìn)行處理，有效避免 了現(xiàn)有方法僅使用標(biāo)量處理方法對(duì)噪聲的敏感性，并且降低了對(duì)結(jié)腸息肉數(shù)據(jù)處理的假陽 性；
[0041] 2、本發(fā)明所提出的方法，通過直接在圖像空間中尋找潛在的等值面，可有效避免 在Marching cubes算法生成的曲面上計(jì)算曲率線所存在的誤差，提高了數(shù)據(jù)處理精度； [0042] 3、本發(fā)明所提出的方法，可自動(dòng)處理圖像數(shù)據(jù)，需要的人機(jī)交互少，提高了智能化 程度；
[0043] 4、本發(fā)明所提出的方法，具有很好的并行性，處理步驟在圖像處理器（圖形處理 器）中并行計(jì)算，適合使用圖像處理器進(jìn)行加速，可有效提高處理速度。
[0044]本發(fā)明的其它優(yōu)點(diǎn)、目標(biāo)和特征將部分通過下面的說明體現(xiàn)，部分還將通過對(duì)本 發(fā)明的研究和實(shí)踐而為本領(lǐng)域的技術(shù)人員所理解。
【附圖說明】
[0045] 圖1為本發(fā)明的面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù)處理方法的系統(tǒng)處理流程 圖；
[0046] 圖2為曲率線追蹤不意圖，圖2(A)為等值面上的曲率線追蹤圖不；圖2(13)為初始點(diǎn) 在等值面上的投影圖示；
[0047]圖3為種子點(diǎn)移動(dòng)步長(zhǎng)的不意圖；
[0048]圖4為曲率線的散播和投影的示意圖；
[0049]圖5為種子點(diǎn)投影距離的示意圖；
[0050]圖6為基于Winding angle的閉合曲率線篩選的示意圖。
【具體實(shí)施方式】
[0051] 下面結(jié)合附圖對(duì)本發(fā)明做進(jìn)一步的詳細(xì)說明，以令本領(lǐng)域技術(shù)人員參照說明書文 字能夠據(jù)以實(shí)施。
[0052] 應(yīng)當(dāng)理解，本發(fā)明所使用的諸如"具有"、"包含"以及"包括"術(shù)語并不配出一個(gè)或 多個(gè)其它元件或其組合的存在或添加。
[0053]本發(fā)明的方法不直接用于疾病診斷，而是用于圖像數(shù)據(jù)的處理過程。
[0054]如圖1-6所示，本發(fā)明提供了一種面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù)處理方法， 包括圖像預(yù)處理、曲率線追蹤、曲率線散播和曲率線特征選擇等4個(gè)主要步驟，其各個(gè)步驟 的計(jì)算過程全部在圖形處理器圖像處理器中進(jìn)行，所有處理步驟在圖像處理器（圖形處理 器）中并行計(jì)算，適合使用圖像處理器進(jìn)行加速，可有效提高處理速度。
[0055] 其中，步驟一、采集結(jié)腸部位的CT三維圖像數(shù)據(jù)f，并將該三維圖像數(shù)據(jù)f讀入到圖 像處理器中，分割出結(jié)腸壁圖像，使用Marching Cubes算法計(jì)算生成分布結(jié)腸壁的所述擬 合等值面曲面，將所述擬合等值面曲面映射到所述三維圖像的坐標(biāo)系中，所述擬合等值面 曲面由若干個(gè)頂點(diǎn)和邊組成，將頂點(diǎn)坐標(biāo)和頂點(diǎn)序號(hào)一一對(duì)應(yīng)，然后進(jìn)行曲率線起點(diǎn)選擇 和篩選，具體操作步驟為：
[0056] (1)、計(jì)算三角面片任意頂點(diǎn)X的主曲率最大主曲率mi和最小主曲率m2。將擬合等值 面曲面劃分為若干區(qū)域，每一個(gè)區(qū)域記作R 1G = I，…，I)。在圖像處理器中設(shè)置I個(gè)block， 記作bi;每一個(gè)block中包含T個(gè)線程r，記作rt(t = 1，2,…，T)。每一個(gè)線程對(duì)應(yīng)的寄存器中 存儲(chǔ)當(dāng)前頂點(diǎn)Vj以及與Vj存在邊連接的頂點(diǎn)坐標(biāo)，這些頂點(diǎn)記作Vx ; Vx、νχ+ι以及當(dāng)前頂點(diǎn)Vj 構(gòu)成了若干個(gè)三角面片，記作
[0057] (11)對(duì)于當(dāng)前頂點(diǎn)Vj，首先，每一個(gè)線程計(jì)算的單位法向量，記作
，并計(jì) 算Vj的法向量
。其中，丐定義為
N是以Vj作為頂點(diǎn)的三角面 片的數(shù)量；
[0058] (12)計(jì)算Vj所在曲線％^的曲_
嗔中，_%&表示頂點(diǎn)Vj和Vx組 成的向1
_示向量Sa的長(zhǎng)度，X是向量乘法。其中，曲線是向量ndPev,$構(gòu)成 的平面與曲面相交所形成的曲線；
[0059] (13)計(jì)算nj和其他向量(經(jīng)過頂點(diǎn)Vj的邊)構(gòu)成平面和曲面所構(gòu)成曲線的曲率。其 經(jīng)過力的所有曲線的曲率構(gòu)成集合&V，取集合中的最大曲率和最小曲率作為該點(diǎn)力的主曲 率，最大曲率記作，最小曲率記作;
[0060] (14)每一個(gè)線程將頂點(diǎn)對(duì)應(yīng)的頂點(diǎn)坐標(biāo)、最大曲率和最小曲率組成一位數(shù)組，放 入寄存器中；
[0061] (2)確定種子點(diǎn)。其具體步驟為：
[0062] (21)每一個(gè)線程判斷所計(jì)算頂點(diǎn)的主曲率，當(dāng)每個(gè)頂點(diǎn)的主曲率滿足以下條件 時(shí)，該頂點(diǎn)作為候選種子點(diǎn)s'：
[0063]
[0064] 其中，mi是最大曲率，m2是最小曲率，mmean是該頂點(diǎn)的平均曲率，mmean = (Ι?1+Π12)/2， mt是預(yù)先選取的閾值，計(jì)算候選種子點(diǎn)s'在其鄰域半徑r范圍內(nèi)所有相鄰頂點(diǎn)Sx的主曲率， 若Sx的最大主曲率和最小主曲率符號(hào)相反，則將候選種子點(diǎn)S'刪除，若S x的最大主曲率和最 小主曲率符號(hào)相同，則將該候選種子點(diǎn)S'作為種子點(diǎn)S放入堆棧S中，假定這些選取體現(xiàn)結(jié) 腸息肉特征的種子點(diǎn)有W個(gè)，其中，線程數(shù)量和候選種子點(diǎn)數(shù)量一致；
[0065] (22)每一個(gè)線程將一個(gè)種子點(diǎn)序號(hào)和坐標(biāo)放入局部存儲(chǔ)器中，并通過顯存將全部 種子點(diǎn)的序號(hào)和坐標(biāo)復(fù)制到內(nèi)存中。
[0066] 步驟二、逐步移動(dòng)所述種子點(diǎn)并將移動(dòng)后的種子點(diǎn)影到圖像真實(shí)等值面上，連接 種子點(diǎn)及其在圖像真實(shí)等值面上的投影點(diǎn)形成W條曲率線;具體包括：
[0067] (I)、計(jì)算每個(gè)種子點(diǎn)s的最大主曲率ki和最小主曲率k2;主曲率ki和k 2的大小為ki = 4"http://^//，k2 = _A2/7/g//。具體計(jì)算過程為:每一個(gè)圖像處理器線程拾取一個(gè)種子點(diǎn)s坐 標(biāo)，計(jì)算種子點(diǎn)s所在圖像灰度曲面的主曲率。其中，g是局部梯度方向，定義為 辦％, 定義Hessian矩陣H為
[0068]
[0069]旋轉(zhuǎn)H，使H的其中一個(gè)坐標(biāo)軸和g重合。旋轉(zhuǎn)后的H記作Hr，Hr可寫為
[0070]
[0071] 其中，（g，u，v)是X的局部坐標(biāo)系。心和\2是二維矩陣Ht的特征值。主曲率1^和1?的方 向?yàn)棣?Ρλ 2所對(duì)應(yīng)的特征向量的方向。其中，可以使用線性插值的方法計(jì)算局部梯度，也可 以使用三次插值方法；
[0072] (II)、在每一個(gè)種子點(diǎn)s對(duì)應(yīng)的真實(shí)等值面的切平面上，圖像處理器線程中的s在 UV平面上沿著最大曲率方向逐步移動(dòng)，其移動(dòng)步長(zhǎng)為L(zhǎng)，其定義為/· = {2?\?(\ + ε, )，其中， lu(i = l，2)是真實(shí)等值面的主曲率，^是設(shè)定的閾值;一般情況下，L = Lmax/10;同時(shí)，我們進(jìn) 一步限定L不得超過0.3毫米;Lmax是三維體素中最長(zhǎng)對(duì)角線的長(zhǎng)度;種子點(diǎn)的在切方向上的 移動(dòng)過程可參見附圖2(A)兒和^的關(guān)系如圖3所示；
[0073] (III)、對(duì)于S中的種子點(diǎn)s，其初始位置為P，該位置處對(duì)應(yīng)的體素灰度值為Ul，由 圖像處理器線程來計(jì)算該位置處的梯度向量#，若Ul小于真實(shí)等值面的體素灰度值U2,則 定義投影向量萬為J，若Ul大于U2，則辦=-愛；
[0074] (VI)、將s在切平面上沿著最大主曲率移動(dòng)后的位置記作P⑴，P⑴沿著D的方向逐 步移動(dòng)，每一步移動(dòng)的步長(zhǎng)為d，其中，d定義為d = Lmax/20，在每一步移動(dòng)之后，重新獲取當(dāng) 前位置Sm的體素灰度，并記錄移動(dòng)前后的體素灰度分別是i前和i后；
[0075] (V)、若（itrU2)X(i后-U2K0,則停止移動(dòng)s，并計(jì)算sm在真實(shí)等值面上的投影位置 P，，
[0076]
[0077] 將P'放入堆棧S中，若種子點(diǎn)在切平面上移動(dòng)的補(bǔ)償LSLmax/lOO,或使用了四階 Runge-Kutta方法追蹤曲率線，則種子點(diǎn)在真實(shí)等值面上的投影過程可以忽略，在曲率線追 蹤停止之前，將投影位置P '作為種子點(diǎn)不斷放入堆棧S中，且將投影位置P '作為種子點(diǎn)回到 步驟(1)，連接在每個(gè)種子點(diǎn)s及其相應(yīng)的若干個(gè)投影位置P'形成一條曲率線；
[0078] (IV)、當(dāng)移動(dòng)后的種子點(diǎn)和其余的曲率線距離小于預(yù)設(shè)的閾值，或者種子點(diǎn)移動(dòng) 次數(shù)達(dá)到了預(yù)設(shè)的上限，則曲率線終止，最終形成W條曲率線，將堆棧S中的種子點(diǎn)信息存儲(chǔ) 到圖像處理器的共享存儲(chǔ)器或顯存中。種子點(diǎn)在圖像等值面的投影過程可參見圖2(B)。
[0079] 步驟三、將每一條所述曲率線進(jìn)行散播，生成曲率線集合。曲率線散播的具體步驟 包括：
[0080] 步驟A、在已經(jīng)生成的曲率線集合中任意選取一條曲率線1，選取1上的一個(gè)種子點(diǎn) s，真實(shí)等值面曲面的切線方向記作T，由g和T定義一個(gè)法平面N;如圖4中的虛線框所示； [0081 ] 步驟B、定義兩個(gè)向量VjPV2 ,V1 = TXg ,V1 = gXT，并開辟兩個(gè)圖像處理器線程"和 t2，tl和t2分別在Vl和V2的方向上分別放置兩個(gè)新增種子點(diǎn)Sl和S2，S1和S2到S的距離為ds，
[0082]
[0083] 若81到已有種子點(diǎn)的距離大于預(yù)設(shè)閾值τ，則將81放入堆棧S;若82到已有種子點(diǎn)的 距離大于預(yù)設(shè)閾值τ，則將 82放入堆棧S;不滿足條件的新增種子點(diǎn)對(duì)應(yīng)的線程終止計(jì)算，其 中，lu(i = l，2)是主曲率，￡1是設(shè)定值;一般情況下，ei = Lmax/10是一個(gè)最優(yōu)選擇；
[0084] 步驟C、分別以81和82為起始點(diǎn)，重復(fù)步驟一和二，生成新的曲率線 終止條件為:移動(dòng)后的種子點(diǎn)和其余的曲率線距離小于預(yù)設(shè)的閾值vt，或者種子點(diǎn)移動(dòng)次 數(shù)達(dá)到了預(yù)設(shè)的上限，其中，Vt為
[0085]
[0086]參量ε2是關(guān)于E1的函數(shù)，可表達(dá)為
[0087]
[0088] 步驟D、曲率線UPl2生成后，刪除線程t#Pt2,對(duì)于UP1 2,在UPl2上分別選取一 個(gè)種子點(diǎn)，并重復(fù)步驟B和C，生成對(duì)應(yīng)曲率線;候選種子點(diǎn)數(shù)量=2X新生成曲率線數(shù)量;使 用二叉樹結(jié)構(gòu)存儲(chǔ)種子點(diǎn)，同一深度的二叉樹節(jié)點(diǎn)中的種子點(diǎn)在圖像處理器中是并行計(jì)算 的；
[0089] 步驟E、重復(fù)上述步驟D，直至所生成的全部候選種子點(diǎn)到已有種子點(diǎn)的距離小于 預(yù)設(shè)閾值τ，曲率線散播結(jié)束，以上步驟全部完成后，生成曲率線集合Ln;然后使用最大曲 率，重復(fù)以上全部步驟，生成曲率線集合L m，Lm是所有沿著最大曲率方向生成的曲率線集 合，Ln是所有沿著最小曲率方向生成的曲率線集合，Ln和LmS所有曲率線的集合。
[0090] 步驟四、在所述曲率線集合中篩選出體現(xiàn)結(jié)腸息肉形狀特征的特征曲率線，并將 所述特征曲率線繪制到所述擬合等值面上，以突顯結(jié)腸息肉。篩選特征曲率線的步驟包括：
[0091] 步驟I、為了有效篩選位于結(jié)腸息肉，構(gòu)造雙曲面百分比參數(shù)HP，W = $ X1 * 其中，nh是位于雙曲面上的曲率線1的種子點(diǎn)，N為曲率線1上的全部種子點(diǎn)，若曲率線1上的 種子點(diǎn)位于某一個(gè)潛在的雙曲面，貝1J其主曲率的符號(hào)相反;在曲率線選擇時(shí)，一般選擇HP最 大的3~5條曲率線。同時(shí)，要求被選擇的曲率線的HP需大于50%;
[0092]步驟J、使用Winding angle表征曲率線的閉合性，記種子點(diǎn)s在移動(dòng)前后分別記作 slPs+，經(jīng)過s的切平面記作T，向量s-s在T上的投影向量記作（si) '，向量SS+在T上的投影 向量記作（Ss+) '，Winding angle是投影向量(si) '到（SS+) '的夾角，如圖6所示，對(duì)篩選后 的每條曲率線計(jì)算其Winding angle，選擇Winding angle最大的若干條曲率線，并且進(jìn)一 步要求曲率線的起點(diǎn)和終點(diǎn)小于預(yù)設(shè)的距離閾值;通過直接在圖像空間中尋找潛在的等值 面，可有效避免在Marching cubes算法生成的曲面上計(jì)算曲率線所存在的誤差，提高了數(shù) 據(jù)處理精度；
[0093] 步驟K、計(jì)算篩選后的曲率線的平均半徑，其定義為曲率線平均中心到曲率線上各 種子點(diǎn)的平均距離，取平均曲率大于預(yù)設(shè)平均曲率閾值rt的曲率線作為特征曲率線，一般 情況下，r t等于5毫米，最后將所述特征曲率線繪制到所述擬合等值面上，以突顯結(jié)腸息肉。
[0094] 由上所述，本發(fā)明提出的方法，使用曲率線這一向量形式對(duì)圖像數(shù)據(jù)進(jìn)行處理，有 效避免了現(xiàn)有方法僅使用標(biāo)量處理方法對(duì)噪聲的敏感性，并且降低了對(duì)結(jié)腸息肉數(shù)據(jù)處理 的假陽性；同時(shí)，通過直接在圖像空間中尋找潛在的等值面，可有效避免在Marching cubes 算法生成的曲面上計(jì)算曲率線所存在的誤差，提高了數(shù)據(jù)處理精度;并且，本發(fā)明的方法可 自動(dòng)處理圖像數(shù)據(jù)，需要的人機(jī)交互少，提高了智能化程度;進(jìn)一步的，本發(fā)明所提出的方 法，具有很好的并行性，處理步驟在圖像處理器（圖形處理器）中并行計(jì)算，適合使用圖像處 理器進(jìn)行加速，可有效提高處理速度。
[0095] 盡管本發(fā)明的實(shí)施方案已公開如上，但其并不僅僅限于說明書和實(shí)施方式中所列 運(yùn)用，它完全可以被適用于各種適合本發(fā)明的領(lǐng)域，對(duì)于熟悉本領(lǐng)域的人員而言，可容易地 實(shí)現(xiàn)另外的修改，因此在不背離權(quán)利要求及等同范圍所限定的一般概念下，本發(fā)明并不限 于特定的細(xì)節(jié)和這里示出與描述的圖例。
【主權(quán)項(xiàng)】
1. 一種面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù)處理方法，其特征在于，在圖形處理器中 進(jìn)行以下步驟： 步驟一、采集結(jié)腸部位的三維圖像數(shù)據(jù)，并對(duì)所述圖像數(shù)據(jù)進(jìn)行圖像預(yù)處理，得到結(jié)腸 壁的擬合等值面曲面，在所述擬合等值面上選取體現(xiàn)結(jié)腸息肉特征的W個(gè)種子點(diǎn)； 步驟二、逐步移動(dòng)所述種子點(diǎn)并將移動(dòng)后的種子點(diǎn)影到圖像真實(shí)等值面上，連接種子 點(diǎn)及其在圖像真實(shí)等值面上的投影點(diǎn)形成W條曲率線； 步驟三、將每一條所述曲率線進(jìn)行散播，生成曲率線集合； 步驟四、在所述曲率線集合中篩選出體現(xiàn)結(jié)腸息肉形狀特征的特征曲率線，并將所述 特征曲率線繪制到所述擬合等值面上，以突顯結(jié)腸息肉。2. 如權(quán)利要求1所述的面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù)處理方法，其特征在于，所 述步驟一中，圖形處理器讀入結(jié)腸部位的CT三維圖像f，分割出結(jié)腸壁圖像，使用Marching Cubes算法計(jì)算生成分布結(jié)腸壁的所述擬合等值面曲面。3. 如權(quán)利要求2所述的面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù)處理方法，其特征在于，所 述步驟一中，將所述擬合等值面曲面映射到所述三維圖像的坐標(biāo)系中，所述擬合等值面曲 面由若干個(gè)頂點(diǎn)和邊組成，將頂點(diǎn)坐標(biāo)和頂點(diǎn)序號(hào) 對(duì)應(yīng)。4. 如權(quán)利要求3所述的面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù)處理方法，其特征在于，所 述步驟一中還包括:計(jì)算每個(gè)頂點(diǎn)的最大主曲率mi和最小主曲率m2，當(dāng)每個(gè)頂點(diǎn)的主曲率滿 足以下條件時(shí)，該頂點(diǎn)作為候選種子點(diǎn)s' ：其中，該頂點(diǎn)的平均曲率，mm_ = (mi+m2) /2，mt是預(yù)先選取的閾值，計(jì)算候選種子 點(diǎn)s'在其鄰域半徑r范圍內(nèi)所有相鄰頂點(diǎn)SX的主曲率，若SX的最大主曲率和最小主曲率符號(hào) 相反，則將候選種子點(diǎn)s'刪除，若s x的最大主曲率和最小主曲率符號(hào)相同，則將該候選種子 點(diǎn)S '作為種子點(diǎn)s放入堆棧S中。5. 如權(quán)利要求4所述的面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù)處理方法，其特征在于，所 述步驟二中包括以下步驟： (I) 、計(jì)算每個(gè)種子點(diǎn)s的最大主曲率ki和最小主曲率k2; (II) 、在每一個(gè)種子點(diǎn)s對(duì)應(yīng)的真實(shí)等值面的切平面上，沿著最大曲率方向逐步移動(dòng)， 其移動(dòng)步長(zhǎng)為L(zhǎng)，其中，ki(i = l，2)是真實(shí)等值面的主曲率，^ 是設(shè)定的閾值；(III) 、對(duì)于S中的種子點(diǎn)s，其初始位置為P，該位置處對(duì)應(yīng)的體素灰度值為U1，由圖像 處理器線程來計(jì)算該位置處的梯度向量:|，若U1小于真實(shí)等值面的體素灰度值U2,則定義 投影向量D為！，若U1大于U 2,則乃=-|; (VI)、將s在切平面上沿著最大主曲率移動(dòng)后的位置記作P⑴，P⑴沿著D的方向逐步移 動(dòng)，每一步移動(dòng)的步長(zhǎng)為d，其中，d定義為d = Lmax/20,在每一步移動(dòng)之后，重新獲取當(dāng)前位 置Sm的體素灰度，并記錄移動(dòng)前后的體素灰度分別是i前和i后； (V)、若（itrU2) X (i后-U2)〈0，則停止移動(dòng)s，并計(jì)算Sm在真實(shí)等值面上的投影位置P'，將P'放入堆棧S中，若種子點(diǎn)在切平面上移動(dòng)的補(bǔ)償LSLmax/100,或使用了四階 Runge-Kutta方法追蹤曲率線，則種子點(diǎn)在真實(shí)等值面上的投影過程可以忽略，在曲率線追 蹤停止之前，將投影位置P '作為種子點(diǎn)不斷放入堆棧S中，且將投影位置P '作為種子點(diǎn)回到 步驟(1)，連接在每個(gè)種子點(diǎn)s及其相應(yīng)的若干個(gè)投影位置P'形成一條曲率線； (IV)、當(dāng)移動(dòng)后的種子點(diǎn)和其余的曲率線距離小于預(yù)設(shè)的閾值，或者種子點(diǎn)移動(dòng)次數(shù) 達(dá)到了預(yù)設(shè)的上限，則曲率線終止，將堆棧S中的種子點(diǎn)信息存儲(chǔ)到圖像處理器的共享存儲(chǔ) 器或顯存中。6. 如權(quán)利要求5所述的面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù)處理方法，其特征在于，所 述步驟三中，曲率線散播的具體步驟包括： 步驟A、在已經(jīng)生成的曲率線集合中任意選取一條曲率線1，選取1上的一個(gè)種子點(diǎn)s，真 實(shí)等值面曲面的切線方向記作T，由g和T定義一個(gè)法平面N; 步驟B、定義兩個(gè)向量Vi和￥2，￥1 = 1\8，'\^ = 8\1'，并開辟兩個(gè)圖像處理器線程1：1和七2,1：1 和t2分別在Vl和V2的方向上分別放置兩個(gè)新增種子點(diǎn)S1和S2，S1和S2到S的距離為ds，若81到已有種子點(diǎn)的距離大于預(yù)設(shè)閾值I，則將81放入堆棧S;gs2到已有種子點(diǎn)的距離 大于預(yù)設(shè)閾值τ，則將82放入堆棧S;不滿足條件的新增種子點(diǎn)對(duì)應(yīng)的線程終止計(jì)算，其中，lu (i = l，2)是主曲率，ει是設(shè)定值； 步驟C、分別以sjPsA起始點(diǎn)，重復(fù)步驟一和二，生成新的曲率線11和12，11和1 2的終止 條件為:移動(dòng)后的種子點(diǎn)和其余的曲率線距離小于預(yù)設(shè)的閾值vt，或者種子點(diǎn)移動(dòng)次數(shù)達(dá) 到了預(yù)設(shè)的上限，其中，vt為參量ε2是關(guān)于^的函數(shù)，可表達(dá)為步驟D、曲率線1#Ρ12生成后，刪除線程t#Pt2,對(duì)于1#Ρ12,在ljPl 2上分別選取一個(gè)種 子點(diǎn)，并重復(fù)步驟B和C，生成對(duì)應(yīng)曲率線;使用二叉樹結(jié)構(gòu)存儲(chǔ)種子點(diǎn)，同一深度的二叉樹 節(jié)點(diǎn)中的種子點(diǎn)在圖像處理器中是并行計(jì)算的； 步驟E、重復(fù)上述步驟D，直至所生成的全部候選種子點(diǎn)到已有種子點(diǎn)的距離小于預(yù)設(shè) 閾值τ，曲率線散播結(jié)束，以上步驟全部完成后，生成曲率線集合Ln;然后使用最大曲率，重 復(fù)以上全部步驟，生成曲率線集合L m，Lm是所有沿著最大曲率方向生成的曲率線集合，匕是 所有沿著最小曲率方向生成的曲率線集合，匕和匕為所有曲率線的集合。7. 如權(quán)利要求6所述的面向虛擬結(jié)腸鏡的結(jié)腸息肉圖像數(shù)據(jù)處理方法，其特征在于，所 述步驟四中，篩選特征曲率線的步驟包括： 步驟I、構(gòu)造雙曲面百分比參數(shù)HP，，其中，nh是位于雙曲面上的曲率線1 的種子點(diǎn)，N為曲率線1上的全部種子點(diǎn)，若曲率線1上的種子點(diǎn)位于某一個(gè)潛在的雙曲面， 則其主曲率的符號(hào)相反； 步驟J、使用Winding angle表征曲率線的閉合性，記種子點(diǎn)s在移動(dòng)前后分別記作?Γ和s +，經(jīng)過s的切平面記作Τ，向量iTs在Τ上的投影向量記作（s^s) '，向量ss+在Τ上的投影向量記 作（ss+) '，Winding angle是投影向量(s^s) '到（ss+) '的夾角，對(duì)篩選后的每條曲率線計(jì)算 其Winding angle，選擇Winding angle最大的若干條曲率線，并且進(jìn)一步要求曲率線的起 點(diǎn)和終點(diǎn)小于預(yù)設(shè)的距離閾值； 步驟K、計(jì)算篩選后的曲率線的平均半徑，其定義為曲率線平均中心到曲率線上各種子 點(diǎn)的平均距離，取平均曲率大于預(yù)設(shè)平均曲率閾值rt的曲率線作為特征曲率線。
【文檔編號(hào)】G06T7/00GK106056583SQ201610347934
【公開日】2016年10月26日
【申請(qǐng)日】2016年5月24日
【發(fā)明人】趙凌霄, 周志勇, 佟寶同, 陳光強(qiáng), 耿辰, 胡冀蘇, 劉燕, 戴亞康
【申請(qǐng)人】中國科學(xué)院蘇州生物醫(yī)學(xué)工程技術(shù)研究所