基于gpu的快速三維ct迭代重建方法
【專利摘要】本發(fā)明公開基于GPU的快速三維CT迭代重建方法，包括以下步驟：步驟一、數(shù)據(jù)輸入，需要將投影數(shù)據(jù)并輸入到數(shù)據(jù)輸入模塊；步驟二、預處理，對投影數(shù)據(jù)作相關的預處理，并將處理后的數(shù)據(jù)以及與重建相關的參數(shù)傳入到GPU；步驟三、為正/反投影，包括正投影步驟、記錄信息步驟和修正步驟、反投影步驟操作，在GPU上實現(xiàn)了正投影步驟和反投影步驟，分別得到正投影系統(tǒng)矩陣和反投影系統(tǒng)矩陣；步驟四、變量更新，按具體迭代算法所對應的迭代公式更新當前的變量值；步驟五、迭代終止，判斷當前迭代是否滿足迭代的終止條件；步驟六、結果輸出，將迭代結果輸出。本發(fā)明只需計算一次射線與體素的相交情況，減少了計算系統(tǒng)矩陣所需的計算量，加快迭代重建的速度。
【專利說明】
基于GPU的快速三維CT迭代重建方法
技術領域
[0001] 本發(fā)明屬于X-射線斷層成像(CT)技術領域，特別涉及基于GPU的快速三維CT迭代 重建方法。
【背景技術】
[0002] 計算機斷層成像技術(Computed Tomography，CT)己經(jīng)廣泛應用于人體組織成像、 工業(yè)無損檢測等領域。CT重建算法可以分為迭代重建算法和解析重建算法兩大類。解析重 建算法的優(yōu)點是算法簡單、重建速度快，缺點是對數(shù)據(jù)的完備性要求較高，并且投影數(shù)據(jù)中 的噪聲對重建圖像的質量影響非常大。相比解析重建算法，迭代重建算法的優(yōu)點是重建圖 像質量高，且適用于各種形式的采集數(shù)據(jù)，即使是在有限角度投影數(shù)據(jù)的情況下，也能重建 出較好的圖像，典型的迭代算法有SART、ISRA、WLS和EM-ML等。但是由于迭代重建算法計算 量大、收斂速度慢、重建時間長，目前已成為制約其廣泛應用的最大瓶頸。
[0003] 迭代重建算法的迭代過程可分解為四個步驟:正投影、修正、反投影和體素更新。 在迭代算法中，與正投影、反投影采用同一驅動方式相比，正投影使用射線驅動方式、反投 影使用體素驅動方式的組合能有效去除環(huán)狀偽影。因此，在迭代算法中，正投影采用基于射 線的驅動方式計算系統(tǒng)矩陣；而反投影采用基于體素驅動方法計算系統(tǒng)矩陣。兩種驅動方 式的區(qū)別在于前者只考慮射線穿過體素的密度，后者還需考慮周圍的體素對當前射線穿過 的體素的影響。在代碼實現(xiàn)過程中，通常將正投影和修正步驟合并，一起進行處理；反投影 和體素更新步驟合并，同時處理。
[0004]在迭代重建算法其中，正投影和反投影的計算量非常大，占迭代重建算法總的計 算量的比例很高，因此提高迭代重建算法的運行速度，關鍵在于提高正、反投影的計算效 率。其中，求解射線與體素的相交情況占用了正投影和反投影的絕大部分的計算量。對于 傳統(tǒng)的方法，正投影和反投影是分別進行的正投影需要確定射線與每個體素的相交情況； 在進行反投影的過程，該操作還得重新進行一次。
[0005] 目前，GPU(圖形處理器，英語:Graphics Processing Unit，縮寫:GPU)技術的出現(xiàn) 減少了傳統(tǒng)迭代算法的運行時間，因為射線驅動方式的正投影和使用體素驅動的反投影是 支持GPU并行加速的。但重建效率依然較低，難以滿足重建效率要求高的應用。因此，提出基 于GPU的快速三維CT迭代重建方法和系統(tǒng)，在理論和實際應用中都具有較高的價值。

【發(fā)明內容】

[0006] 本發(fā)明提出了基于GPU的快速三維CT迭代重建方法，主要是針對傳統(tǒng)三維CT迭代 重建所需計算時間長的問題，提高重建效率。
[0007] 本發(fā)明通過以下技術方案實現(xiàn)：
[0008] 基于GPU的快速三維CT迭代重建方法，包括以下步驟：
[0009] 步驟一、數(shù)據(jù)輸入，需要通過CT設備采集投影數(shù)據(jù)并輸入到數(shù)據(jù)輸入模塊，以及獲 取該投影數(shù)據(jù)時的具體指模體位置和尺寸、旋轉中心的位置、射線源到旋轉中心的距離、探 測器尺寸、探測器到旋轉中心的距離和投影角度；
[0010] 步驟二、預處理，預處理模塊對投影數(shù)據(jù)作相關的預處理，并將處理后的數(shù)據(jù)以及 與重建相關的參數(shù)傳入到GPU;
[0011] 步驟三、為正/反投影，正/反投影模塊進行包括正投影步驟、記錄信息步驟和修正 步驟、反投影步驟操作，該模塊在GPU上實現(xiàn)了正投影步驟和反投影步驟，分別得到正投影 系統(tǒng)矩陣和反投影系統(tǒng)矩陣，具體包括:首先，選出一部分待計算正投影的射線，通過多線 程分別對這些選中的射線進行正投影，記錄每條射線正投影的結果，以及射線與體素相交 的信息;然后，在選中射線的正投影步驟和信息記錄步驟完成之后，根據(jù)記錄的信息進行修 正步驟和反投影步驟操作，其中，修正過程只與投影數(shù)據(jù)和每條射線的正投影有關，根據(jù) 之前正投影步驟記錄的信息可以得到當前迭代過程中每條射線的投影值，從而可以快速地 計算每條射線的修正后投影值，在反投影過程中，根據(jù)之前正投影步驟記錄的信息可以快 速地得到每條射線穿過體素的位置和長度，從而可以根據(jù)每條射線與體素的相交情況，快 速地將每條射線的修正后投影值分配到相交位置處相鄰的8個體素上，最后，繼續(xù)從剩余的 射線中選取一部分射線重復之前的操作，直到所有的射線都被計算過而終止循環(huán)；
[0012] 步驟四、變量更新，變量更新模塊按具體迭代算法所對應的迭代公式更新當前的 變量值；
[0013] 步驟五、迭代終止，迭代終止模塊判斷當前迭代是否滿足迭代的終止條件，如果滿 足迭代的終止條件，終止迭代，輸出迭代結果到結果輸出模塊，否則，繼續(xù)之前的正/反投影 操作和變量更新操作；
[0014] 步驟六、結果輸出，結果輸出模塊將迭代結果輸出，包括存儲和顯示迭代重建后的 數(shù)據(jù)。
[0015] 進一步，所述的反投影步驟可以通過多線程操作進行加速，每個線程都需要分配 一塊數(shù)據(jù)來存儲各體素的值。
[0016] 進一步，所述的正/反投影步驟中的修正步驟由所使用的具體迭代重建算法決定。
[0017] 進一步，所述的迭代重建算法為EM-ML算法，以下是EM-ML算法的迭代公式：
[0019]其中，模體內有N個體素，在第k次迭代中，每個體素位置處的劑量值為<，利用M條 射線照射模體，投影數(shù)據(jù)為y，在公式中，;為正投影過程
為修正因子， 對的值進行修正，將修正后的值進行反投影
，最后利用反投影值更新當前 的xf值，在正投影過程中，使用基于射線的驅動方式計算系統(tǒng)矩陣a^，而反投影過程采用基 于體素的驅動方式計算系統(tǒng)矩陣，兩種驅動方式的差別在于:射線j穿過體素 i的線長度 為1，對于基于射線的驅動方式，值為1和體素 i的密度值的乘積;對于基于體素的驅動方 式，bij值為1和1中間位置處對應密度值的乘積，通過對中間位置處周圍8個體素的密度值進 行插值得到該密度值，在正/反投影過程中使用非匹配的系統(tǒng)矩陣是為了減少重建圖像中 的環(huán)狀偽影，先計算正投影;，記錄每條射線的正投影值和的值及其位置;然后， 計算每條射線的修正因子，在前面兩個過程中，可以用多線同時計算多條射線的正投影值 及修正因子;最后，利用插值算法，將修正因子與的乘積賦給體素 i及周圍其它7個體素， 在反投影過程中同樣可以使用多線程進行加速。
[0020] 本發(fā)明與現(xiàn)有技術相比，優(yōu)點如下：
[0021] 1、正/反投影過程中，只需計算一次射線與體素的相交情況，減少了計算系統(tǒng)矩陣 所需的計算量，加快迭代重建的速度。
[0022] 2、正/反投影模塊將正投影步驟和反投影步驟合并，將減少確定射線與體素相交 情況的次數(shù)，加快迭代重建的速度。
【附圖說明】：
[0023]下面結合附圖對本發(fā)明的【具體實施方式】作進一步詳細說明。
[0024]圖1是本發(fā)明基于GPU的快速三維CT迭代重建方法的具體實施流程圖；
[0025]圖2是本發(fā)明正/反投影的具體實施示意圖；
[0026]圖3是本發(fā)明正投影和記錄信息的具體實施示意圖；
[0027]圖4是本發(fā)明修正和反投影的具體實施示意圖；
[0028]圖5是本發(fā)明基于GPU的快速三維CT迭代重建系統(tǒng)框圖。
【具體實施方式】
[0029] 為使本發(fā)明實施的目的、特點和優(yōu)點能夠更為明顯易懂，下面結合附圖對本發(fā)明 的具體實施例作詳細的說明。
[0030] 本發(fā)明提供了基于GPU的快速三維CT迭代重建方法和系統(tǒng)，如圖5所示，基于GPU的 快速三維CT迭代重建系統(tǒng)100包括數(shù)據(jù)輸入模塊1，預處理模塊2,正/反投影模塊3,變量更 新模塊4,迭代終止模塊5,結果輸出模塊6;如圖1~圖5所示，為該系統(tǒng)的具體實施流程圖：
[0031] 步驟S101為數(shù)據(jù)輸入，需要通過CT設備采集投影數(shù)據(jù)并輸入到數(shù)據(jù)輸入模塊1，以 及獲取該投影數(shù)據(jù)時的模體位置和尺寸、旋轉中心的位置、射線源到旋轉中心的距離、探測 器尺寸、探測器到旋轉中心的距離和投影角度；
[0032]步驟S102為預處理，預處理模塊2對投影數(shù)據(jù)作相關的預處理，例如，為了降低投 影數(shù)據(jù)中的噪聲，可以對投影數(shù)據(jù)進行非線性濾波、基于貝葉斯統(tǒng)計理論的濾波;為了克服 重建后圖像中的條形偽影、環(huán)狀偽影、金屬偽影和杯形偽影，可以采用基于字典學習、形態(tài) 分量分析等新興的圖像處理方法對投影數(shù)據(jù)進行處理，并將處理后的數(shù)據(jù)以及與重建相關 的參數(shù)傳入到GPU，所述處理后的數(shù)據(jù)以及與重建相關的參數(shù)具體指模體位置和尺寸、旋轉 中心的位置、射線源到旋轉中心的距離、探測器尺寸、探測器到旋轉中心的距離和投影角 度，由于這些處理后的數(shù)據(jù)和參數(shù)是在CPU上，這些數(shù)據(jù)還需傳入到GPU上才能使用GPU進行 加速；
[0033] 步驟S103為正/反投影，正/反投影模塊3進行包括正投影步驟S301、記錄信息步驟 S302和修正步驟S401、反投影步驟S402操作（如圖2所示），該模塊在GPU上實現(xiàn)正投影步驟 S301和反投影步驟S301，分別計算得到正投影系統(tǒng)矩陣和反投影系統(tǒng)矩陣，如圖2、圖3、圖4 所示，該正/反投影模塊3內進行了正投影步驟S301和反投影步驟S401，分別得到正投影系 統(tǒng)矩陣和反投影系統(tǒng)矩陣，具體包括：
[0034]首先，選出一部分待計算正投影的射線，通過多線程分別對這些選中的射線進行 正投影(如圖3的S301)，記錄每條射線正投影的結果，以及射線與體素相交的信息(如圖3的 S302)，
[0035] 然后，在選中射線的正投影步驟S301和信息記錄步驟S302完成之后，根據(jù)記錄的 信息進行修正步驟S401和反投影步驟S402操作，
[0036]其中，修正過程只與投影數(shù)據(jù)和每條射線的正投影有關，根據(jù)之前正投影步驟記 錄的信息可以得到當前迭代過程中每條射線的投影值，從而可以快速地計算每條射線的修 正后投影值，在反投影過程中，根據(jù)之前正投影步驟記錄的信息可以快速地得到每射線穿 過體素的位置和長度，從而可以根據(jù)每條射線與體素的相交情況，快速地將每條射線的修 正后投影值分配到相交位置處相鄰的8個體素上，這時的反投影操作避免了重新確定射線 和體素的相交情況，減少了算法的計算量，反投影操作可以通過多線程操作進行加速，但是 同時進行基于體素的操作容易產(chǎn)生"寫沖突"，每個線程都需要分配一塊數(shù)據(jù)來存儲各體素 的值，因此反投影操作不能啟用過多的線程，
[0037]最后，繼續(xù)從剩余的射線中選取一部分射線重復之前的操作，直到所有的射線都 被計算過而終止循環(huán)；
[0038]所述的正/反投影模塊3中正/反投影過程的修正步驟S401由所使用的具體迭代重 建算法決定，為了明顯易懂，以EM-ML算法為例進行詳細說明，以下是EM-ML算法的迭代公 式：
[0040] 其中，模體內有N個體素，在第k次迭代中，每個體素位置處的劑量值為，利用M條 射線照射模體，投影數(shù)據(jù)為y，在公式中，;為正投影過程：
為修正因子， 對的值進行修正，將修正后的值進行反投影
最后利用反投影值更新當前 的<值，在正投影過程中，使用基于射線的驅動方式計算系統(tǒng)矩陣a^，而反投影過程采用基 于體素的驅動方式計算系統(tǒng)矩陣，兩種驅動方式的差別在于:射線j穿過體素 i的線長度 為1，對于基于射線的驅動方式，值為1和體素 i的密度值的乘積;對于基于體素的驅動方 式，bij值為1和1中間位置處對應密度值的乘積，通過對中間位置處周圍8個體素的密度值進 行插值得到該密度值，在正/反投影過程中使用非匹配的系統(tǒng)矩陣是為了減少重建圖像中 的環(huán)狀偽影，先計算正投影，記錄每條射線的正投影值和的值及其位置;然后， 計算每條射線的修正因子，在前面兩個過程中，可以用多線同時計算多條射線的正投影值 及修正因子(如圖3所示）；最后，利用插值算法，將修正因子與的乘積賦給體素 i及周圍其 它7個體素，在反投影過程中同樣可以使用多線程進行加速；
[0041] 步驟S104為變量更新，變量更新模塊4按具體迭代算法所對應的迭代公式更新當 前的;cf值；
[0042]步驟S105為迭代終止，迭代終止模塊5判斷當前迭代是否滿足迭代的終止條件，如 果滿足迭代的終止條件，終止迭代，輸出迭代結果到結果輸出模塊，否則，繼續(xù)之前的正/反 投影操作和變量更新操作；
[0043]步驟S106為結果輸出，結果輸出模塊6將迭代結果輸出，包括存儲和顯示迭代重建 后的數(shù)據(jù)。
[0044] 本發(fā)明書中未作詳細描述的內容屬于本領域專業(yè)技術人員公知的現(xiàn)有技術。
[0045] 以上所述僅是本發(fā)明的優(yōu)選實施方式，應當指出，對于本技術領域的普通技術人 員來說，在不脫離本發(fā)明原理的前提下，還可以做出若干的改進和潤飾，這些改進和潤飾 也應視為本發(fā)明的保護范圍。
【主權項】
1. 基于GPU的快速三維CT迭代重建方法，其特征在于:包括以下步驟： 步驟一、數(shù)據(jù)輸入，需要通過CT設備采集投影數(shù)據(jù)并輸入到數(shù)據(jù)輸入模塊，以及獲取該 投影數(shù)據(jù)時的具體指模體位置和尺寸、旋轉中心的位置、射線源到旋轉中心的距離、探測器 尺寸、探測器到旋轉中心的距離和投影角度； 步驟二、預處理，預處理模塊對投影數(shù)據(jù)作相關的預處理，并將處理后的數(shù)據(jù)以及與重 建相關的參數(shù)傳入到GPU; 步驟三、為正/反投影，正/反投影模塊進行包括正投影步驟、記錄信息步驟和修正步 驟、反投影步驟操作，該模塊在GPU上實現(xiàn)了正投影步驟和反投影步驟，分別得到正投影系 統(tǒng)矩陣和反投影系統(tǒng)矩陣，具體包括:首先，選出一部分待計算正投影的射線，通過多線程 分別對這些選中的射線進行正投影，記錄每條射線正投影的結果，以及射線與體素相交的 信息;然后，在選中射線的正投影步驟和信息記錄步驟完成之后，根據(jù)記錄的信息進行修正 步驟和反投影步驟操作，其中，修正過程只與投影數(shù)據(jù)和每條射線的正投影有關，根據(jù)之前 正投影步驟記錄的信息可以得到當前迭代過程中每條射線的投影值，從而可以快速地計算 每條射線的修正后投影值，在反投影過程中，根據(jù)之前正投影步驟記錄的信息可以快速地 得到每條射線穿過體素的位置和長度，從而可以根據(jù)每條射線與體素的相交情況，快速地 將每條射線的修正后投影值分配到相交位置處相鄰的8個體素上，最后，繼續(xù)從剩余的射線 中選取一部分射線重復之前的操作，直到所有的射線都被計算過而終止循環(huán)； 步驟四、變量更新，變量更新模塊按具體迭代算法所對應的迭代公式更新當前的變量 值； 步驟五、迭代終止，迭代終止模塊判斷當前迭代是否滿足迭代的終止條件，如果滿足迭 代的終止條件，終止迭代，輸出迭代結果到結果輸出模塊，否則，繼續(xù)之前的正/反投影操 作和變量更新操作； 步驟六、結果輸出，結果輸出模塊將迭代結果輸出，包括存儲和顯示迭代重建后的數(shù) 據(jù)。2. 如權利要求1所述的基于GPU的快速三維CT迭代重建方法，其特征在于:所述的反投 影步驟可以通過多線程操作進行加速，每個線程都需要分配一塊數(shù)據(jù)來存儲各體素的值。3. 如權利要求1所述的基于GPU的快速三維CT迭代重建方法，其特征在于:所述的正/反 投影步驟中的修正步驟由所使用的具體迭代重建算法決定。4. 如權利要求3所述的基于GPU的快速三維CT迭代重建方法，其特征在于:所述的迭代 重建算法為EM-ML算法，以下是EM-ML算法的迭代公式：其中，模體內有N個體素，在第k次迭代中，每個體素位置處的劑量值為xf，利用M條射線 照射模體，投影數(shù)據(jù)為y，在公式中，為正投影過程，！為修正因子，對 的值進行修正，將修正后的值進行反投影>最后利用反投影值更新當前的xf 值，在正投影過程中，使用基于射線的驅動方式計算系統(tǒng)矩陣a^，而反投影過程采用基于體 素的驅動方式計算系統(tǒng)矩陣by，兩種驅動方式的差別在于:射線j穿過體素 i的線長度為1， 對于基于射線的驅動方式，值為1和體素 i的密度值的乘積;對于基于體素的驅動方式， bij值為1和1中間位置處對應密度值的乘積，通過對中間位置處周圍8個體素的密度值進行 插值得到該密度值，在正/反投影過程中使用非匹配的系統(tǒng)矩陣是為了減少重建圖像中的 環(huán)狀偽影，先計算正投影，記錄每條射線的正投影值和的值及其位置;然后，計 算每條射線的修正因子，在前面兩個過程中，可以用多線同時計算多條射線的正投影值及 修正因子;最后，利用插值算法，將修正因子與的乘積賦給體素 i及周圍其它7個體素，在 反投影過程中同樣可以使用多線程進行加速。
【文檔編號】G06T17/00GK105931280SQ201610217047
【公開日】2016年9月7日
【申請日】2016年3月29日
【發(fā)明人】張權, 張鵬程, 劉祎, 桂志國, 韓躍平, 楊民
【申請人】中北大學