專利名稱：一種旋轉(zhuǎn)x射線造影圖像迭代重建方法
技術(shù)領(lǐng)域：
本發(fā)明屬于計算機技術(shù)領(lǐng)域，涉及ー種旋轉(zhuǎn)X射線造影圖像迭代重建方法。
背景技術(shù)：
旋轉(zhuǎn)X射線冠狀動脈造影成像技術(shù)是繼雙平面血管造影成像以來又ー項得到廣泛關(guān)注的冠狀動脈成像技木。通過該技木，醫(yī)師能夠從多個角度提供完整，準確地觀察血管的形狀以及運動方式，對腦血管腫瘤，冠心病的診療具有重要意義。但僅憑旋轉(zhuǎn)圖像序列，難以對冠狀動脈的具體結(jié)構(gòu)有一個直觀準確的認識，因此血管造影的三維重建成為目前學(xué)術(shù)領(lǐng)域的研究熱點，也是醫(yī)療器械廠商亟待解決的重要問題。旋轉(zhuǎn)X射線造影技術(shù)的成像幾何同三維重建同錐束CT類似。由于旋轉(zhuǎn)造影成像過程中伴隨著心臟運動，為了獲得準確的血管三維結(jié)構(gòu)，需要采用迭代重建算法。但是由于該成像系統(tǒng)的空間分辨率高于一般的 錐束CT，該重建問題的計算量遠遠大于錐束CT圖象重建。迭代算法的運算主要集中在投影、反投影兩個環(huán)節(jié)，投影、反投影的快速計算問題成為目前旋轉(zhuǎn)X射線血管造影成像三維重建問題的主要瓶頸，投影矩陣的使用可以大大加速投影，反投影操作的運算速度，但受到硬件條件限制，目前的計算機內(nèi)存存儲能力不足以存儲旋轉(zhuǎn)造影系統(tǒng)投影矩陣。

發(fā)明內(nèi)容
技術(shù)問題本發(fā)明提供了一種對硬件系統(tǒng)要求較低，能夠大幅度提高投影、反投影速度，便于并性計算優(yōu)化，可快速準確獲得血管三維結(jié)構(gòu)的旋轉(zhuǎn)X射線造影圖像迭代重建方法。技術(shù)方案本發(fā)明的旋轉(zhuǎn)X射線造影圖像迭代重建方法，包括以下步驟I)從旋轉(zhuǎn)X射線造影設(shè)備讀取掃描數(shù)據(jù)文件，保存投影序列圖像并記錄如下參數(shù)射線源到探測板距離SDD，射線源到C臂旋轉(zhuǎn)中心距離S0D，各投影采樣的旋轉(zhuǎn)角度，投影圖像像素邊長h，投影圖像長度U個像素単位和寬度V個像素単位，所述投影圖像尺寸與ニ維投影空間尺寸一致；確定ニ維投影空間坐標(biāo)軸U、V，所述坐標(biāo)軸U、V分別平行于ニ維投影空間的長、寬方向，根據(jù)用戶的精度需求，分別設(shè)定三維圖像空間長度為X個體素単位、寬度為Y個體素単位、高度為Z個體素単位，三維圖像空間體素邊長1，確定三維圖像空間坐標(biāo)軸x、y、z，所述坐標(biāo)軸x、y、z分別平行于三維圖像空間的長、寬、高方向，且均通過三維圖像空間的中心位置；2)根據(jù)步驟I)中設(shè)定的三維圖像空間長寬高對三維圖像空間進行降采樣操作，具體方法為選定降采樣倍數(shù)dif，所述降采樣倍數(shù)dif為能被三維圖像空間長度X、寬度Y、高度Z分別整除的整數(shù)，將三維圖像空間中由dif*dif*diff體素組成的、邊長為dif*l的立方體內(nèi)的體素，歸并為降采樣體素，對整個三維圖像空間完成歸并處理后，按照歸并順序排列所述降采樣體素，得到降采樣三維圖像空間，降采樣體素邊長為三維圖像空間體素邊長的dif倍;根據(jù)步驟I)中記錄的投影圖像長寬，對ニ維投影空間進行降采樣操作，具體方法為選定降采樣倍數(shù)dpf，所述降采樣倍數(shù)為能被ニ維投影空間長度U和寬度V分別整除的整數(shù)，將ニ維投影空間中由dpf*dpf個像素組成的、邊長為dpf*h的正方形內(nèi)的像素，歸并為投影空間降采樣像素，對整個ニ維投影空間完成歸并處理后，按照歸并順序排列所述投影空間降采樣像素，得到降采樣ニ維投影空間，投影空間降采樣像素邊長為ニ維投影空間像素邊長的Clpf倍；根據(jù)步驟I)中記錄的投影圖像長寬，對步驟I)獲得的投影序列圖像進行降采樣操作，得到降采樣序列投影圖像，具體方法為對投影序列中每個投影圖像，將其中由cU*dpf個像素組成的、邊長為dpf*h的正方形內(nèi)的像素，歸并為投影圖像降采樣像素，并且對所述正方形內(nèi)的像素值進行累加求和，作為投影圖像降采樣像素值，對整個投影圖像完成歸并處理后，按照歸并順序排列所述投影圖像降采樣像素，得到降采樣投影圖像，投影圖像降采樣像素邊長為投影圖像像素邊長的dpf倍；3)根據(jù)步驟I)中所記錄的探測板距離SDD，射線源到C臂旋轉(zhuǎn)中心距離S0D，初始采樣旋轉(zhuǎn)角度，針對步驟2)得到的降采樣三維圖像空間和降采樣ニ維投影空間，利用距離驅(qū)動算法，構(gòu)造初始掃描方向低分辨率體素索引投影矩陣，并將得到的初始掃描方向低分 辨率體素索引投影矩陣以三元組存放方式保存于內(nèi)存設(shè)備中；4)針對步驟2)提供的降采樣三維圖像空間，構(gòu)造各旋轉(zhuǎn)角度下降采樣三維圖像空間旋轉(zhuǎn)矩陣R，具體方法為首先根據(jù)圖像旋轉(zhuǎn)角度位置關(guān)系，利用線性插值算法構(gòu)造并存儲降采樣圖像空間切面層旋轉(zhuǎn)矩陣R°，隨后根據(jù)降采樣三維圖像空間旋轉(zhuǎn)矩陣R的分塊對
稱性，利用公式R =得到降采樣圖像空間旋轉(zhuǎn)矩陣R ；5)對步驟2)得到的降采樣投影序列圖像采用頂帽濾波方法進行濾波處理，對濾波處理結(jié)果ニ值化得到第一分割結(jié)果，對步驟2)得到的降采樣投影序列圖像采用Frangi血管濾波方法進行濾波處理，對濾波處理結(jié)果ニ值化得到第二分割結(jié)果，對第一分割結(jié)果和第二分割結(jié)果求并集，得到低分辨率投影序列圖像血管分割結(jié)果；6)利用步驟3)所得到的三元組格式存儲的初始掃描方向低分辨率體素索引投影矩陣和步驟4)所得到的各旋轉(zhuǎn)角度降采樣圖像空間旋轉(zhuǎn)矩陣，對步驟5)所得到的低分辨率投影序列圖像血管分割結(jié)果進行反投影操作具體方法為對低分辨率投影圖像血管分割結(jié)果，根據(jù)像素編號排列拉伸低分辨率投影圖像血管分割結(jié)果，得到低分辨率投影圖像血管分割結(jié)果向量，將所得向量同初始掃描方向低分辨率體素索引投影矩陣的轉(zhuǎn)置矩陣相乘得到中間低分辨率反投影結(jié)果；將該中間低分辨率反投影結(jié)果同對應(yīng)旋轉(zhuǎn)角度降采樣圖像空間旋轉(zhuǎn)矩陣的轉(zhuǎn)置矩陣相乘，得到該旋轉(zhuǎn)角度下的低分辨率反投影結(jié)果；7)利用步驟6)所得到的低分辨率反投影結(jié)果確定低分辨率三維圖像空間血管掩摸，具體方法為當(dāng)投影數(shù)目小于等于5吋，對各旋轉(zhuǎn)角度下低分辨率反投影結(jié)果求交集，作為低分辨率三維圖像空間血管掩模；當(dāng)投影數(shù)目大于5時，利用對各旋轉(zhuǎn)角度下低分辨率反投影結(jié)果求和并進行閾值劃分，作為低分辨率三維圖像空間血管掩模；8)對步驟7)得到的低分辨率三維圖像空間血管掩模，根據(jù)步驟2)所設(shè)定的降采樣倍數(shù)dif進行升采樣操作，得到三維圖像空間血管掩模，具體方法為將低分辨率三維圖像空間血管掩模的每ー個低分辨率體素在x,y,z方向等分為dif份，拆分成為dif*dif*dif個高分辨率體素，將該低分辨率體素的值作為其拆分成的高分辨率體素的值，并按照拆分順序排列高分辨率體素，當(dāng)?shù)头直媛嗜S圖像空間血管掩模的所有低分辨率體素完成上述操作后，即得到三維圖像空間血管掩模；9)根據(jù)步驟I)中所記錄的射線源到探測板距離SDD，射線源到C臂旋轉(zhuǎn)中心距離S0D，初始采樣旋轉(zhuǎn)角度，投影圖像長度U、寬度V，投影圖像像素邊長h，設(shè)定的三維圖像空間長度X、寬度Y、高度Z，三維圖像空間體素邊長1，以及步驟8)所得的到三維圖像空間血管掩模，利用距離驅(qū)動算法計算各旋轉(zhuǎn)角度下的完整投影矩陣；10)根據(jù)步驟I)所記錄的投影序列圖像，根據(jù)步驟9)所得到的各旋轉(zhuǎn)角度下的完整投影矩陣，完成三維血管結(jié)構(gòu)重建，具體方法為101)設(shè)定ー個重建結(jié)果向量，所述重建結(jié)果向量長度為圖像空間體素總個數(shù)XXYXZ，重建結(jié)果向量元素值全部為O;設(shè)定ー個步長向量，步長向量長度為圖像空間體素總個數(shù)XXYXZ，步長向量元素值全部為O ;設(shè)定ー個單元向量，単元向量長度為圖像空間體素總個數(shù)XXYXZ，単元向量元素值全部為I ; 102)對于每個投影角度，進行如下操作將單元向量與當(dāng)前投影角度下的完整投影矩陣相乘，將相乘結(jié)果同當(dāng)前投影角度下的完整投影矩陣的轉(zhuǎn)置矩陣相乘，得到當(dāng)前投影角度的單方向步長向量，將該單方向步長向量累加至步長向量；對所有投影角度完成上述操作后，將累加完成的步長向量記錄保存為迭代步長向量；103)執(zhí)行重建結(jié)果向量迭代更新步驟100至300次，完成指定次數(shù)的重建結(jié)果向量迭代更新步驟后，將重建結(jié)果向量保存為重建結(jié)果三維血管體數(shù)據(jù)，所述重建結(jié)果向量迭代更新的方法如下對于每個投影角度，進行如下操作將當(dāng)前的重建結(jié)果向量同當(dāng)前投影角度下的完整投影矩陣相乘，得到當(dāng)前方向臨時投影向量，用當(dāng)前方向投影圖像向量減去所述當(dāng)前方向臨時投影向量，得到當(dāng)前方向臨時投影差向量，將所述當(dāng)前方向臨時投影差向量同當(dāng)前投影角度下的完整投影矩陣的轉(zhuǎn)置矩陣相乘，得到當(dāng)前方向差異反投影向量，用所述當(dāng)前方向差異反投影向量點除以迭代步長向量，得到當(dāng)前方向迭代步進向量，將當(dāng)前方向迭代步進向量累加至重建結(jié)果向量；對所有角度完成上述操作后，完成一次重建結(jié)果向量更新。本發(fā)明中，步驟3)中利用距離驅(qū)動算法，構(gòu)造初始掃描方向低分辨率體素索引投影矩陣，并將得到的初始掃描方向低分辨率體素索引投影矩陣以三元組存放方式保存于內(nèi)存設(shè)備中的具體方法為31)設(shè)定發(fā)射源的空間坐標(biāo)為(_sSQD，0，0)，設(shè)定探測板中心點坐標(biāo)為(sSDD-sSQD，0，O), Ssod表示射線源到C臂旋轉(zhuǎn)中心距離，Ssdd表示射線源到探測板距離；32)對所有降采樣體素作如下操作從步驟31)設(shè)定的發(fā)射源向序號為jd的體素各角點引出射線，各射線將與降采樣ニ維投影空間相交，其八個交點的連線所包圍區(qū)域面積記<_，所述包圍區(qū)域同序號為id的像素相交面積為Sd' μ，則au=sd' JiZsdj-即為低分辨率投影矩陣第id行，第jd列的元素值；33)將所得到的所有具有非零值的元素的行號、列號以及元素值保存，得到三元組格式存儲的初始掃描方向低分辨率體素索引投影矩陣。本發(fā)明的步驟4)中構(gòu)造并存儲降采樣圖像空間切面層旋轉(zhuǎn)矩陣R°的具體方法為記投影采樣的旋轉(zhuǎn)角度為α，將三維圖像空間按照ζ方向體素単位分組，共計Z個切面層，選取任一切面層，將其中的每個體素圍繞Z軸旋轉(zhuǎn)-a角度，記錄旋轉(zhuǎn)后的位置的鄰域體素編號Itl, I1, 12，I3，按照線性插值方法確定矩陣第k行，第Itl列，第I1列，第I2列，第I3列的四個元素的元素值，其中k為體素的序號，將所得到的所有具有非零值的元素的行號，列號以及元素值保存，得到三元組格式存儲的降采樣圖像空間切面層旋轉(zhuǎn)矩陣R°。本發(fā)明中，步驟9)中利用距離驅(qū)動算法計算各旋轉(zhuǎn)角度下的完整投影矩陣的具體過程為91)設(shè)定發(fā)射源的空間坐標(biāo)為(-SsqdXcoS a，-sS0DXsina , O)，設(shè)定探測板中心點坐豐不為((Ssdd-Ssod) X cos Ct，(sSDD—sS0D) X sin ct，O);92)對所有三維圖像空間體素作如下操作如果在步驟8所得到的三維圖像空間血管掩模中，序號為j的體素沒有非零值，則不做任何操作直接進入下一個體素的操作，如果具有非零值，則根據(jù)步驟91)設(shè)定的發(fā)射源向該體素各角點引出射線，各射線將與降采樣 ニ維投影空間相交，其八個交點的連線所包圍區(qū)域面積記為ち，同序號為i的像素相交面積為t, ノも即為投影矩陣第i行，第j列的元素值；93)將所得到的所有具有非零值的元素的行號，列號以及元素值保存，得到三元組格式存儲的α角度方向投影矩陣。本發(fā)明是一種旋轉(zhuǎn)X射線造影圖像迭代重建方法，該方法分為3個階段第一階段構(gòu)造低分辨率投影矩陣，并將完整矩陣拆解為單一角度矩陣和旋轉(zhuǎn)矩陣2個分量進行簡化存儲。在該階段首先在単一投影方向上，利用體素索引投影算法構(gòu)造低分辨率投影矩陣；之后構(gòu)造3D體數(shù)據(jù)旋轉(zhuǎn)矩陣，利用該旋轉(zhuǎn)矩陣自身対稱性進行簡化存儲。第二階段在第一階段得到低分辨率投影矩陣的基礎(chǔ)上，進行進一步基于投影內(nèi)容的簡化。在該階段，首先對采集投影數(shù)據(jù)進行降采樣及快速分割；之后利用低分辨率投影矩陣進行反投影，確定3D熱區(qū)，其后對該區(qū)域進行升采樣。最后利用升采樣區(qū)域構(gòu)造投影矩陣掩摸，結(jié)合系統(tǒng)參數(shù)計算精確投影矩陣。第三階段在第二階段得到的精確投影矩陣的基礎(chǔ)上進行旋轉(zhuǎn)X射線造影系統(tǒng)的三維血管重建。在該階段，首先利用単元向量計算步長向量；之后通過一定次數(shù)的迭代步驟，反復(fù)計算修正當(dāng)前結(jié)果向量同投影數(shù)據(jù)之間的差異，井根據(jù)步長向量設(shè)定每次的具體修正量，在迭代步驟完成后獲得三維血管重建數(shù)據(jù)。有益效果本發(fā)明同現(xiàn)有技術(shù)方法相比，具有如下優(yōu)點本方法確定的投影矩陣構(gòu)造方法大幅降低了矩陣存儲空間，一方面，利用稀疏矩陣存儲技術(shù)，可以降低所需存儲空間，另一方面，同一般稀疏矩陣存儲方式相比，考慮了血管掩模后，所需存儲空間可進一歩大幅降低。以如下數(shù)據(jù)為例ニ維投影大小為512*512，三維圖像大小為256*256*256 :單ー投影角度完整投影矩陣大小為512*512*256*256*256，以每個元素值占用4字節(jié)存儲空間為例，共需要存儲空間16384GB ;如果使用稀疏矩陣存儲，由于平均每個體素對應(yīng)相關(guān)像素為12個左右，則矩陣非零元素個數(shù)為256*256*256*9，每個非零元素的存儲需要耗費4字節(jié)元素值，以及4字節(jié)行列號，共需要存儲空間約為I. 7GB ;利用本發(fā)明提出的算法，在考慮圖像掩模的情況下，通常血管掩模的非零元素個數(shù)在2*105左右，每個非零元素對應(yīng)12個相關(guān)像素，共需要存儲空間約為20Mb字節(jié)。同一般稀疏矩陣相比，本方法所需存儲空間降至1%左右，利用獲得的投影矩陣同投影向量相乘，則完成投影操作，利用獲得的投影矩陣的轉(zhuǎn)置，同圖像向量相乘，則完成反投影操作。同一般稀疏矩陣相比本方法確定的投影矩陣構(gòu)造方法大幅提高了投影、反投影計算速度，記錄通過矩陣相乘運算完成投影、反投影操作所需時間，理論上，矩陣尺寸減小至1%，則速度提高應(yīng)為100倍，實際應(yīng)用當(dāng)中，由于數(shù)據(jù)傳遞，數(shù)據(jù)交互等底層原因，實驗證明速度提升約為20-40倍。獲得簡化后的投影矩陣后，可以利用其完成旋轉(zhuǎn)X射線造影系統(tǒng)的三維血管重建，得到準確的三維血管結(jié)構(gòu)。


圖I是本發(fā)明方法的旋轉(zhuǎn)造影系統(tǒng)投影生成示意圖。圖2a是本發(fā)明方法的距離驅(qū)動像素索引投影矩陣幾何示意圖。圖2b是本發(fā)明方法的距離驅(qū)動像素索引投影矩陣元素權(quán)重計算示意圖。 圖3是本發(fā)明方法的ニ值化投影分割反投影疊加示意圖。圖4是本發(fā)明方法的整體流程圖。圖5是本發(fā)明方法步驟3)的子流程圖。圖6是本發(fā)明方法步驟4)中構(gòu)造并存儲降采樣圖像空間切面層旋轉(zhuǎn)矩陣R°的子流程圖。圖7是本發(fā)明方法步驟9)中利用距離驅(qū)動算法計算各旋轉(zhuǎn)角度下的完整投影矩陣的子流程圖。
具體實施例方式本發(fā)明的旋轉(zhuǎn)X射線造影圖像迭代重建方法，包括下列步驟I)從旋轉(zhuǎn)X射線造影設(shè)備讀取掃描數(shù)據(jù)文件，保存投影序列圖像并記錄如下參數(shù)射線源到探測板距離SDD，射線源到C臂旋轉(zhuǎn)中心距離S0D，各投影采樣的旋轉(zhuǎn)角度，投影圖像像素邊長h，投影圖像像長度U個像素単位和寬度V個像素単位，投影圖像尺寸與ニ維投影空間像素尺寸一致，每個象素單位為長度為h的正方形，確定ニ維投影空間坐標(biāo)軸U，V，坐標(biāo)軸U，V分別平行于ニ維投影空間的長，寬方向，根據(jù)用戶的精度需求，分別設(shè)定三維圖像空間長度為X個體素単位、寬度為Y個體素単位、高度為Z個體素単位，每個體素單位為邊長I的立方體，確定三維圖像空間坐標(biāo)軸X，y, z，坐標(biāo)軸X，y, z分別平行于三維圖像空間的長、寬、高方向，且均通過三維圖像空間的中心位置，也即坐標(biāo)原點為三維圖像空間的中心位置。系統(tǒng)模擬圖如圖I所示，圖中S為發(fā)射源位置，O為旋轉(zhuǎn)中心，D為探測板中心位置。射線源到探測板距離=SD，射線源到C臂旋轉(zhuǎn)中心距離=S0，記圖像空間體素尺寸IX 1X1，圖像空間包含體素個數(shù)XXYXZ ;記投影圖像像素尺寸為hXh，投影空間包含像素個數(shù)UX V;記向量X/)和X軸夾角為a，每次以O(shè)為中心旋轉(zhuǎn)角度為Λ a，總投影角度數(shù)為K。2)根據(jù)步驟I)中設(shè)定的三維圖像空間長寬高對三維圖像空間進行降采樣操作，具體方法為選定降采樣倍數(shù)dif，所述降采樣倍數(shù)dif為能被三維圖像空間長度X、寬度Y、高度Z分別整除的整數(shù)，通?？梢赃x擇2，3，4，將三維圖像空間中由dif*dif*dif個體素組成的、邊長為dif*l的立方體內(nèi)的體素，歸并為降采樣體素，對整個三維圖像空間完成歸并處理后，按照歸并順序排列所述降采樣體素，得到降采樣三維圖像空間，降采樣體素邊長為三維圖像空間體素邊長的dif倍，記Xif = X/dif, Yif = Y/dif, Zif = Z/dif則降采樣圖像空間體素尺寸IdifX IdifXldif,降采樣圖像空間包含體素個數(shù)為Xif*Yif*Zif ；根據(jù)步驟I)中記錄的投影圖像長寬，對ニ維投影空間進行降采樣操作，具體方法為選定降采樣倍數(shù)dpf，所述降采樣倍數(shù)為能被ニ維投影空間長度U和寬度V分別整除的整數(shù)，通常可以選擇2，3 ，4，將ニ維投影空間中由dpf*dpf個像素組成的、邊長為dpf*h的正方形內(nèi)的像素，歸并為投影空間降采樣像素，對整個ニ維投影空間完成歸并處理后，按照歸并順序排列所述投影空間降采樣像素，得到降采樣ニ維投影空間，投影空間降采樣像素邊長為ニ維投影空間像素邊長的Clpf倍，記Upf = U/dpf, Vpf = V/dpf則降采樣投影空間像素尺寸hdpf Xhdpf,降采樣投影空間包含像素個數(shù)UpfXVpf ；根據(jù)步驟I)中記錄的投影圖像長寬，對步驟I)獲得的投影序列圖像進行降采樣操作，得到降采樣序列投影圖像，具體方法為對投影序列中每個投影圖像，將其中由cU*dpf個像素組成的、邊長為dpf*h的正方形內(nèi)的像素，歸并為投影圖像降采樣像素，并且對所述正方形內(nèi)的像素值進行累加求和，作為投影圖像降采樣像素值，對整個投影圖像完成歸并處理后，按照歸并順序排列所述投影圖像降采樣像素，得到降采樣投影圖像，投影圖像降采樣像素邊長為投影圖像像素邊長的Clpf倍，降采樣投影圖像像素尺寸hdpf Xhdpf,降采樣投影圖像包含像素個數(shù)Upf X Vpf。3)根據(jù)步驟I)中所記錄的探測板距離SDD，射線源到C臂旋轉(zhuǎn)中心距離S0D，初始采樣旋轉(zhuǎn)角度，針對步驟2)得到的降采樣三維圖像空間和降采樣ニ維投影空間，利用距離驅(qū)動算法，構(gòu)造初始掃描方向低分辨率體素索引投影矩陣，具體方法為取初始方向a=0，設(shè)定發(fā)射源的空間坐標(biāo)為(_sS(I)，0，0)，探測板中心坐標(biāo)為(sSDD-sS(I)，0，0)，按照X，y, z的順序遍歷降采樣圖像空間體素。如圖2 (a)所示，對于序號為(VX，Vy，VZ)的體素，記其編號jd = vzXXifXYif+vyXXif+vx,作點源S和編號jd體素中心點連線,其延長線與探測板相交點在探測平面的坐標(biāo)為(pud，pvd)，作點源S和編號jd體素各角點連線，其延長線與探測板相交點所包圍矩形面積記為cC. ,.，其面積為<_，若序號為(m，η)的投影像素點，其編號記為id = pvdXU+pud,與,相交，且相交面積為sd/ u，，如圖2 (b)所示,則記低分辨率體素索引投影矩陣矩陣第id行，jd列的元素值~=パ^./Sdj ；將所有具有非零值的元素的行號，列號以及元素值以三元組存放方式保存于內(nèi)存設(shè)備中，得到的初始掃描方向低分辨率體素索引投影矩陣；4)針對步驟2)提供的降采樣三維圖像空間，構(gòu)造各旋轉(zhuǎn)角度下降采樣三維圖像空間旋轉(zhuǎn)矩陣R具體方法為首先對圖像空間沿z軸(旋轉(zhuǎn)軸)進行劃分，根據(jù)圖像旋轉(zhuǎn)角度位置關(guān)系，利用線性插值算法構(gòu)造降采樣圖像空間切面層旋轉(zhuǎn)矩陣R°，該矩陣大小為，XifYifXXifYif為稀疏矩陣，對于所選切面層，可看作ニ維圖像處理，以圖像中心為原點，對于
像素(row，col)，其2D坐標(biāo)分別為(a°，b°)，其中％ ==し"/-記其序
Z ο
號為j = col XXif+row。對應(yīng)稀疏旋轉(zhuǎn)矩陣R°的第j列。根據(jù)2D旋轉(zhuǎn)坐標(biāo)變換公式，經(jīng)過旋轉(zhuǎn)_a角度后該像素點的坐標(biāo)位置為
權(quán)利要求
1.一種旋轉(zhuǎn)X射線造影圖像迭代重建方法，其特征在于，包括以下步驟 1)從旋轉(zhuǎn)X射線造影設(shè)備讀取掃描數(shù)據(jù)文件，保存投影序列圖像并記錄如下參數(shù)射線源到探測板距離SDD，射線源到C臂旋轉(zhuǎn)中心距離SOD，各投影采樣的旋轉(zhuǎn)角度，投影圖像像素邊長h，投影圖像長度U個像素單位和寬度V個像素單位，所述投影圖像尺寸與二維投影空間尺寸一致；確定二維投影空間坐標(biāo)軸U、V，所述坐標(biāo)軸U、V分別平行于二維投影空間的長、寬方向，根據(jù)用戶的精度需求，分別設(shè)定三維圖像空間長度為X個體素單位、寬度為Y個體素單位、高度為Z個體素單位，三維圖像空間體素邊長1，確定三維圖像空間坐標(biāo)軸χ、y、z，所述坐標(biāo)軸x、y、z分別平行于三維圖像空間的長、寬、高方向，且均通過三維圖像空間的中心位置； 2)根據(jù)步驟I)中設(shè)定的三維圖像空間長寬高對三維圖像空間進行降采樣操作，具體方法為選定降采樣倍數(shù)dif，所述降采樣倍數(shù)dif為能被三維圖像空間長度X、寬度Y、高度Z分 別整除的整數(shù)，將三維圖像空間中由dif*dif*dif個體素組成的、邊長為dif*l的立方體內(nèi)的體素，歸并為降采樣體素，對整個三維圖像空間完成歸并處理后,按照歸并順序排列所述降采樣體素，得到降采樣三維圖像空間，降采樣體素邊長為三維圖像空間體素邊長的dif倍； 根據(jù)步驟I)中記錄的投影圖像長寬，對二維投影空間進行降采樣操作，具體方法為選定降采樣倍數(shù)Clpf，所述降采樣倍數(shù)為能被二維投影空間長度U和寬度V分別整除的整數(shù)，將二維投影空間中由dpf*dpf個像素組成的、邊長為dpf*h的正方形內(nèi)的像素，歸并為投影空間降采樣像素，對整個二維投影空間完成歸并處理后，按照歸并順序排列所述投影空間降采樣像素，得到降采樣二維投影空間，投影空間降采樣像素邊長為二維投影空間像素邊長的Clpf倍； 根據(jù)步驟I)中記錄的投影圖像長寬，對步驟I)獲得的投影序列圖像進行降采樣操作，得到降采樣序列投影圖像，具體方法為對投影序列中每個投影圖像，將其中由dpf*dpf個像素組成的、邊長為dpf*h的正方形內(nèi)的像素，歸并為投影圖像降采樣像素，并且對所述正方形內(nèi)的像素值進行累加求和，作為投影圖像降采樣像素值，對整個投影圖像完成歸并處理后，按照歸并順序排列所述投影圖像降采樣像素，得到降采樣投影圖像，投影圖像降采樣像素邊長為投影圖像像素邊長的Clpf倍； 3)根據(jù)步驟I)中所記錄的探測板距離SDD，射線源到C臂旋轉(zhuǎn)中心距離SOD，初始采樣旋轉(zhuǎn)角度，針對步驟2)得到的降采樣三維圖像空間和降采樣二維投影空間，利用距離驅(qū)動算法，構(gòu)造初始掃描方向低分辨率體素索引投影矩陣，并將得到的初始掃描方向低分辨率體素索引投影矩陣以三元組存放方式保存于內(nèi)存設(shè)備中； 4)針對步驟2)提供的降采樣三維圖像空間，構(gòu)造各旋轉(zhuǎn)角度下降采樣三維圖像空間旋轉(zhuǎn)矩陣R，具體方法為首先根據(jù)圖像旋轉(zhuǎn)角度位置關(guān)系，利用線性插值算法構(gòu)造并存儲降采樣圖像空間切面層旋轉(zhuǎn)矩陣R°，隨后根據(jù)降采樣三維圖像空間旋轉(zhuǎn)矩陣R的分塊對稱性，利用公式Λ =到降采樣圖像空間旋轉(zhuǎn)矩陣R ； Z 5)對步驟2)得到的降采樣投影序列圖像采用頂帽濾波方法進行濾波處理，對濾波處理結(jié)果二值化得到第一分割結(jié)果，對步驟2)得到的降采樣投影序列圖像采用Frangi血管濾波方法進行濾波處理，對濾波處理結(jié)果二值化得到第二分割結(jié)果，對第一分割結(jié)果和第二分割結(jié)果求并集，得到低分辨率投影序列圖像血管分割結(jié)果；6)利用步驟3)所得到的三元組格式存儲的初始掃描方向低分辨率體素索引投影矩陣和步驟4)所得到的各旋轉(zhuǎn)角度降采樣圖像空間旋轉(zhuǎn)矩陣，對步驟5)所得到的低分辨率投影序列圖像血管分割結(jié)果進行反投影操作具體方法為對低分辨率投影圖像血管分割結(jié)果，根據(jù)像素編號排列拉伸低分辨率投影圖像血管分割結(jié)果，得到低分辨率投影圖像血管分割結(jié)果向量，將所得向量同初始掃描方向低分辨率體素索引投影矩陣的轉(zhuǎn)置矩陣相乘得到中間低分辨率反投影結(jié)果；將該中間低分辨率反投影結(jié)果同對應(yīng)旋轉(zhuǎn)角度降采樣圖像空間旋轉(zhuǎn)矩陣的轉(zhuǎn)置矩陣相乘，得到該旋轉(zhuǎn)角度下的低分辨率反投影結(jié)果； 7)利用步驟6)所得到的低分辨率反投影結(jié)果確定低分辨率三維圖像空間血管掩模，具體方法為當(dāng)投影數(shù)目小于等于5時，對各旋轉(zhuǎn)角度下低分辨率反投影結(jié)果求交集，作為低分辨率三維圖像空間血管掩模；當(dāng)投影數(shù)目大于5時，利用對各旋轉(zhuǎn)角度下低分辨率反投影結(jié)果求和并進行閾值劃分，作為低分辨率三維圖像空間血管掩模； 8)對步驟7)得到的低分辨率三維圖像空間血管掩模，根據(jù)步驟2)所設(shè)定的降采樣倍數(shù)dif進行升采樣操作，得到三維圖像空間血管掩模，具體方法為將低分辨率三維圖像空 間血管掩模的每一個低分辨率體素在X，1，Z方向等分為dif份，拆分成為dif*dif*dif個高分辨率體素，將該低分辨率體素的值作為其拆分成的高分辨率體素的值，并按照拆分順序排列高分辨率體素，當(dāng)?shù)头直媛嗜S圖像空間血管掩模的所有低分辨率體素完成上述操作后，即得到三維圖像空間血管掩模； 9)根據(jù)步驟I)中所記錄的射線源到探測板距離SDD，射線源到C臂旋轉(zhuǎn)中心距離SOD，初始采樣旋轉(zhuǎn)角度，投影圖像長度U、寬度V，投影圖像像素邊長h，設(shè)定的三維圖像空間長度X、寬度Y、高度Z，三維圖像空間體素邊長1，以及步驟8)所得的到三維圖像空間血管掩模，利用距離驅(qū)動算法計算各旋轉(zhuǎn)角度下的完整投影矩陣； 10)根據(jù)步驟I)所記錄的投影序列圖像，根據(jù)步驟9)所得到的各旋轉(zhuǎn)角度下的完整投影矩陣，完成三維血管結(jié)構(gòu)重建，具體方法為 101)設(shè)定一個重建結(jié)果向量，所述重建結(jié)果向量長度為圖像空間體素總個數(shù)XXYXZ，重建結(jié)果向量元素值全部為O ;設(shè)定一個步長向量，所述步長向量長度為圖像空間體素總個數(shù)XXYXZ，步長向量元素值全部為O ;設(shè)定一個單元向量，所述單元向量長度為圖像空間體素總個數(shù)XXYXZ，單元向量元素值全部為I ; 102)對于每個投影角度，進行如下操作將單元向量與當(dāng)前投影角度下的完整投影矩陣相乘，將相乘結(jié)果同當(dāng)前投影角度下的完整投影矩陣的轉(zhuǎn)置矩陣相乘，得到當(dāng)前投影角度的單方向步長向量，將該單方向步長向量累加至步長向量；對所有投影角度完成上述操作后，將累加完成的步長向量記錄保存為迭代步長向量； 103)執(zhí)行重建結(jié)果向量迭代更新步驟100至300次，完成指定次數(shù)的重建結(jié)果向量迭代更新步驟后，將重建結(jié)果向量保存為重建結(jié)果三維血管體數(shù)據(jù)，所述重建結(jié)果向量迭代更新的方法如下 對于每個投影角度，進行如下操作將當(dāng)前的重建結(jié)果向量同當(dāng)前投影角度下的完整投影矩陣相乘，得到當(dāng)前方向臨時投影向量，用當(dāng)前方向投影圖像向量減去所述當(dāng)前方向臨時投影向量，得到當(dāng)前方向臨時投影差向量，將所述當(dāng)前方向臨時投影差向量同當(dāng)前投影角度下的完整投影矩陣的轉(zhuǎn)置矩陣相乘，得到當(dāng)前方向差異反投影向量，用所述當(dāng)前方向差異反投影向量點除以迭代步長向量，得到當(dāng)前方向迭代步進向量，將當(dāng)前方向迭代步進向量累加至重建結(jié)果向量；對所有角度完成上述操作后，完成一次重建結(jié)果向量更新。
2.根據(jù)權(quán)利要求I所述的旋轉(zhuǎn)X射線造影圖像迭代重建方法，其特征在于，所述步驟3)中利用距離驅(qū)動算法，構(gòu)造初始掃描方向低分辨率體素索引投影矩陣，并將得到的初始掃描方向低分辨率體素索引投影矩陣以三元組存放方式保存于內(nèi)存設(shè)備中的具體方法為 31)設(shè)定發(fā)射源的空間坐標(biāo)為(_sSQD，0，0)，設(shè)定探測板中心點坐標(biāo)為(sSDD-sSQD，0，0)，Ssod表示射線源到C臂旋轉(zhuǎn)中心距離，Ssdd表示射線源到探測板距離； 32)對所有降采樣體素作如下操作從步驟31)設(shè)定的發(fā)射源向序號為jd的體素各角點引出射線，各射線將與降采樣二維投影空間相交，其八個交點的連線所包圍區(qū)域面積記<，所述包圍區(qū)域同序號為id的像素相交面積為Sd' μ則au=sd' JiZsdj-即為低分辨率投影矩陣第0行，第jd列的元素值；· 33)將所得到的所有具有非零值的元素的行號、列號以及元素值保存，得到三元組格式存儲的初始掃描方向低分辨率體素索引投影矩陣。
3.根據(jù)權(quán)利要求I所述的旋轉(zhuǎn)X射線造影圖像迭代重建方法，其特征在于，所述的步驟4)中構(gòu)造并存儲降采樣圖像空間切面層旋轉(zhuǎn)矩陣R°的具體方法為記投影采樣的旋轉(zhuǎn)角度為α，將三維圖像空間按照ζ方向體素單位分組，共計Z個切面層，選取任一切面層，將其中的每個體素圍繞ζ軸旋轉(zhuǎn)-α角度，記錄旋轉(zhuǎn)后的位置的鄰域體素編號Itl, I1, I2, I3，按照線性插值方法確定矩陣第k行，第Itl列，第I1列，第I2列，第I3列的四個元素的元素值，其中k為體素的序號，將所得到的所有具有非零值的元素的行號，列號以及元素值保存，得到三元組格式存儲的降采樣圖像空間切面層旋轉(zhuǎn)矩陣R°。
4.根據(jù)權(quán)利要求I所述的旋轉(zhuǎn)X射線造影圖像迭代重建方法，其特征在于，所述步驟9)中利用距離驅(qū)動算法計算各旋轉(zhuǎn)角度下的完整投影矩陣的具體過程為 91)設(shè)定發(fā)射源的空間坐標(biāo)為(-SsodXc0Sa，-sS0DXsina, O)，設(shè)定探測板中心點坐標(biāo)為((sSDD-sS0D) Xcos a , (sSDD-sS0D) Xsina , O)； 92)對所有三維圖像空間體素作如下操作如果在步驟8所得到的三維圖像空間血管掩模中，序號為j的體素沒有非零值，則不做任何操作直接進入下一個體素的操作，如果具有非零值，則根據(jù)步驟91)設(shè)定的發(fā)射源向該體素各角點引出射線，各射線將與降采樣二維投影空間相交，其八個交點的連線所包圍區(qū)域面積記為\_，同序號為i的像素相交面積為t/ ji，則au=t' /%即為投影矩陣第i行，第j列的元素值； 93)將所得到的所有具有非零值的元素的行號，列號以及元素值保存，得到三元組格式存儲的α角度方向投影矩陣。
全文摘要
本發(fā)明公開了一種旋轉(zhuǎn)X射線造影圖像迭代重建方法，首先在第一階段構(gòu)造低分辨率投影矩陣，并將完整矩陣拆解為單一角度矩陣和旋轉(zhuǎn)矩陣2個分量進行簡化存儲，然后第二階段在第一階段得到低分辨率投影矩陣的基礎(chǔ)上，進行進一步基于投影內(nèi)容的簡化，最后在第三階段進行三維血管重建本方法采用的經(jīng)過掩模簡化的投影矩陣解決了旋轉(zhuǎn)X射線造影系統(tǒng)投影，反投影計算量過大，計算時間過長的問題，能夠有效獲得三維血管結(jié)構(gòu)，幫助臨床醫(yī)師進行診斷。
文檔編號G06T11/00GK102842141SQ20121022816
公開日2012年12月26日 申請日期2012年7月3日 優(yōu)先權(quán)日2012年7月3日
發(fā)明者胡軼寧, 謝理哲, 沈傲東, 羅立民 申請人:東南大學(xué)