基于cuda架構(gòu)的gpu并行加速雙能譜ct重建方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及X射線CT技術(shù)領(lǐng)域雙能CT解析重建算法，尤其涉及一種基于CUDA架 構(gòu)的GPU并行加速雙能譜CT解析重建方法。
【背景技術(shù)】
[0002] 對(duì)于CT加速領(lǐng)域，采用GPU加速已經(jīng)成為重要的方法，而具體的加速方法是提高 CT產(chǎn)品的處理速度的關(guān)鍵環(huán)節(jié)，既是本領(lǐng)域技術(shù)人員研宄的核心技術(shù)，也直接決定產(chǎn)品的 市場(chǎng)競(jìng)爭(zhēng)力。
[0003] 與單能CT重建不同，雙能CT重建并不是前者的簡(jiǎn)單重復(fù)執(zhí)行，而具有其特殊性。 雙能CT解析重建算法分為后處理重建和預(yù)處理重建兩種算法。其中，后處理重建和單能 CT解析重建方法一致，是對(duì)投影數(shù)據(jù)直接重建，獲得被檢物的高低能線衰減系數(shù)信息，重建 的形狀精度較好，但由于不能從根本上消除射術(shù)硬化效應(yīng)，重建的被檢物材料精度難以保 證；而預(yù)處理重建與后處理算法相比，預(yù)處理重建算法要對(duì)雙能投影數(shù)據(jù)進(jìn)行分解，預(yù)處理 重建算法分為"基效應(yīng)分解模型"和"基材料分解模型"兩種模型，兩種模型是統(tǒng)一的，對(duì)分 解后的基材料（或基效應(yīng)）投影數(shù)據(jù)進(jìn)行重建，再根據(jù)基材料的重建圖像計(jì)算等效原子序 數(shù)和電子密度圖像。這種預(yù)處理算法理論上可以得到不受能譜影響的被檢物的等效原子序 數(shù)和電子密度信息，重建的材料精度較高。雙能CT重建想要同時(shí)獲得上述四種重建圖像信 息，時(shí)間因素是算法能否實(shí)際應(yīng)用的關(guān)鍵。
[0004] 如現(xiàn)有技術(shù)中，公開號(hào)為：101283913,名稱：CT圖像重建的GPU加速方法，該發(fā)明 專利申請(qǐng)公開了一種基于GPU硬件加速的CT圖像重建方法，該發(fā)明軟件部分通常包括基 于GPU的CT數(shù)據(jù)預(yù)處理模塊、基于GPU的CT數(shù)據(jù)濾波模塊、基于GPU的CT正投影模塊和 基于GPU的CT圖像重建反投影模塊，該發(fā)明利用GPU硬件實(shí)現(xiàn)了對(duì)CT圖像重建算法的加 速，重建部分和數(shù)據(jù)處理部分均在GPU上完成，通過對(duì)大數(shù)據(jù)處理采用分塊處理方法，對(duì)每 塊待重建體重建時(shí)只需要用到部分投影數(shù)據(jù)，減少了數(shù)據(jù)傳輸，從而可提高整個(gè)體重建速 度。
[0005] 又如現(xiàn)有技術(shù)中，公開號(hào)為：101596113,名稱：一種CT并行重建系統(tǒng)及成像方法， 該發(fā)明專利申請(qǐng)公開了一種CT并行重建系統(tǒng)及其成像方法，該發(fā)明所述系統(tǒng)包括前端采 樣器、中心節(jié)點(diǎn)、以及與所述中心節(jié)點(diǎn)連接的多個(gè)子節(jié)點(diǎn)，每個(gè)所述子節(jié)點(diǎn)均安裝有圖形 處理器。該發(fā)明所述成像方法包括如下步驟：1)將重建圖像劃分為多個(gè)BLOCK區(qū)域；2)中 心節(jié)點(diǎn)接收前端采樣器所采集的各角度下的原始投影數(shù)據(jù)；3)中心節(jié)點(diǎn)為各節(jié)點(diǎn)分配計(jì) 算任務(wù)，每個(gè)分配了計(jì)算任務(wù)的子節(jié)點(diǎn)計(jì)算一個(gè)或多個(gè)BLOCK區(qū)域的重建值；4)各子節(jié)點(diǎn) 完成所分配的重建計(jì)算；5)中心節(jié)點(diǎn)組合出完整的重建圖像。該發(fā)明采取劃分重建圖像子 區(qū)域的方法充分挖掘CT掃描及重建的并行特性，在數(shù)據(jù)采集的同時(shí)進(jìn)行GPU的并行重建， 重建時(shí)間由現(xiàn)有技術(shù)的分鐘級(jí)降低為毫秒級(jí)，達(dá)到準(zhǔn)實(shí)時(shí)的被測(cè)對(duì)象的斷層圖像重建顯 不〇
[0006] 由此可知，現(xiàn)有技術(shù)基本上都是利用GPU并行實(shí)現(xiàn)單一的CT圖像重建過程，以達(dá) 到加速圖像重建的目的，由于CT重建本身計(jì)算量大，耗時(shí)高，算法速度一直是CT相關(guān)算法 能否在實(shí)際中使用的關(guān)鍵影響因素，對(duì)于雙能CT解析重建算法中預(yù)處理重建算法的特殊 性，現(xiàn)有技術(shù)中并沒有給出有針對(duì)性的完善的并行加速解決方法，基于上述現(xiàn)有技術(shù)中的 不足，這也成為雙能CT解析重建算法發(fā)展過程中亟待解決的問題。

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

[0007] 為了彌補(bǔ)現(xiàn)有技術(shù)的缺陷，本發(fā)明提供一種基于CUDA架構(gòu)的GPU并行加速雙能譜 CT解析重建方法，結(jié)合雙能CT重建的特點(diǎn)，提供一種基于CUDA架構(gòu)利用GPU的并行能力， 進(jìn)行雙能CT快速解析重建同時(shí)獲得四幅重建圖像的方法，S卩：后處理重建的高低能線衰減 系數(shù)圖像和預(yù)處理重建的等效原子序數(shù)和電子密度圖像。
[0008] 為了實(shí)現(xiàn)上述目的，本發(fā)明所述采取的技術(shù)解決方案是：本發(fā)明所述一種基于 CUDA架構(gòu)的GPU并行加速雙能譜CT解析重建方法，所述方法包括以下基本步驟：
[0009] S1 :高低能采樣數(shù)據(jù)快速配準(zhǔn)；
[0010] S2 :并行計(jì)算原始投影數(shù)據(jù)及投影分解和投影加權(quán)；
[0011] S3 :雙能投影或者分解投影同步快速卷積濾波；
[0012] S4 :同時(shí)并行反投影重建幾類圖像；
[0013] 進(jìn)一步的，也可以包括：
[0014] S5 :由分解重建圖像計(jì)算等效原子序數(shù)圖像和電子密度圖像。
[0015] 根據(jù)不同應(yīng)用，上述步驟也可以調(diào)整為S1對(duì)原始投影數(shù)據(jù)進(jìn)行快速配準(zhǔn)；S2雙能 投影分解和投影加權(quán)。
[0016] 進(jìn)一步的，其中所述步驟S1的具體特征為：
[0017] S11 :所述高低能采樣數(shù)據(jù)由雙能譜CT系統(tǒng)采樣獲得，其雙能譜為通過快速切換 射線源產(chǎn)生的真實(shí)雙能譜或者通過使用雙層夾心探測(cè)器獲得的偽雙能譜；
[0018] S12 :所述雙能譜CT系統(tǒng)的掃描模式為圓周掃描或者螺旋掃描；雙能均采用順時(shí) 針掃描或者逆時(shí)針掃描；
[0019] S13 :所述雙能譜CT系統(tǒng)的探測(cè)器為單排探測(cè)器或者多排探測(cè)器或者面陣探測(cè) 器；
[0020] S14 :所述的CT系統(tǒng)相關(guān)參數(shù)在GPU設(shè)備端定義為常量參數(shù)，并在初始化時(shí)賦值， 以便加快GPU內(nèi)核函數(shù)對(duì)其訪問的速度；
[0021] S15:所述的配準(zhǔn)是為保證雙能CT預(yù)處理投影分解重建算法的分解精度，先對(duì)雙 能采樣數(shù)據(jù)進(jìn)行視角配準(zhǔn)，其特征在于：
[0022] S16:所述的雙能快速配準(zhǔn)使用計(jì)算最大互相關(guān)系數(shù)方法，為加速將計(jì)算雙能二維 采樣數(shù)據(jù)的互相關(guān)系數(shù)簡(jiǎn)化為對(duì)采樣數(shù)據(jù)進(jìn)行一維投影，再計(jì)算雙能一維投影的互相關(guān)系 數(shù)；
[0023] S17 :采用GPU并行方式對(duì)高低能采樣進(jìn)行局部探測(cè)器單元數(shù)據(jù)二值化，并對(duì)二值 化后的數(shù)據(jù)進(jìn)行行方向的一維投影；
[0024] S18 :通過雙能一維投影的卷積實(shí)現(xiàn)計(jì)算互相關(guān)系數(shù)加速，使用CUFFT對(duì)雙能采樣 的一維投影變換到頻域進(jìn)行，獲得互相關(guān)系數(shù)數(shù)組，尋找數(shù)組中最大互相關(guān)系數(shù)所在位置 即為配準(zhǔn)需要調(diào)整的行數(shù)。
[0025] 進(jìn)一步地，其中所述步驟S2的具體特征在于：
[0026] S21 :計(jì)算原始雙能投影為加快讀取數(shù)據(jù)的速度，將采樣數(shù)據(jù)和空采明場(chǎng)數(shù)據(jù)都存 儲(chǔ)在紋理存儲(chǔ)器中；
[0027] S22:基于查找表實(shí)現(xiàn)雙能投影分解，初始化時(shí)將預(yù)先離線生成的"投影分解查找 表"存儲(chǔ)到顯存的二維紋理存儲(chǔ)器以便提高讀取速度，該查找表紋理采用非歸一化的浮點(diǎn) 型拾取坐標(biāo)，將高低能投影值作為橫縱坐標(biāo)查找對(duì)應(yīng)的分解投影值，并利用紋理存儲(chǔ)器硬 件實(shí)現(xiàn)的線性模式濾波功能對(duì)讀取的浮點(diǎn)型返回值進(jìn)行插值，實(shí)現(xiàn)快速高精度雙能投影分 解；
[0028] S23:在同一內(nèi)核函數(shù)中按照像素級(jí)并行先計(jì)算高低能原始投影數(shù)據(jù)，再利用高低 能原始投影數(shù)據(jù)到查找表紋理中查找分解投影數(shù)據(jù)，并對(duì)原始投影和分解投影的同時(shí)進(jìn)行 加權(quán)操作。根據(jù)具體應(yīng)用場(chǎng)景不同，加權(quán)值可以離線計(jì)算后，初始化時(shí)導(dǎo)入顯存紋理中，在 加權(quán)操作時(shí)讀取，或者在該內(nèi)核函數(shù)中實(shí)現(xiàn)同步計(jì)算該像素點(diǎn)的加權(quán)值。在同一內(nèi)核函數(shù) 中同時(shí)實(shí)現(xiàn)上述計(jì)算的優(yōu)點(diǎn)在于節(jié)省存取中間結(jié)果的時(shí)間，加快處理速度。
[0029] 進(jìn)一步地，其中所述步驟S3的具體特征在于：
[0030] S31 :濾波過程中使用CUDA架構(gòu)提供的基于GPU的快速傅里葉變換函數(shù)庫(kù)CUFFT;
[0031] S32 :初始化時(shí)將濾波核用FFT變換到頻域，并存儲(chǔ)到一維紋理存儲(chǔ)器中，以便在 頻域運(yùn)算時(shí)加快讀取速度；
[0032] S33 :同步進(jìn)行雙能和雙分解后投影的FFT變換，利用復(fù)數(shù)的實(shí)部和虛部同時(shí)進(jìn)行 高低能投影的FFT變換或者分解后的雙投影的FFT變換；
[0033] S34 :利用FFT同時(shí)處理一批一維離散傅里葉變換，并將變換后的復(fù)數(shù)值存儲(chǔ)到復(fù) 數(shù)紋理以便計(jì)算頻域點(diǎn)積時(shí)加快讀取速度；
[0034] S35:使用同一內(nèi)核函數(shù)同步對(duì)雙能和雙分解后頻域投影進(jìn)行與濾波核的點(diǎn)積計(jì) 算，將濾波后的結(jié)果數(shù)據(jù)綁定到紋理存儲(chǔ)器，以便反投影重建時(shí)加快讀取速度。
[0035] 進(jìn)一步地，其中所述步驟S4的具體特征在于：
[0036] S41 :采用像素驅(qū)動(dòng)的反投影重建，不存在對(duì)存儲(chǔ)重建圖像的全局顯存的"寫沖 突"；
[0037] S42 :同一內(nèi)核函數(shù)中按照像素級(jí)并行進(jìn)行幾類圖像反投影重建，同一線程完成四 種類型重建圖像的同一像素位置的反投影計(jì)算，即同時(shí)完成高低能線衰減系數(shù)圖像和兩種 分解圖像中各自的同一像素位置的重建；
[0038] S43 :重建圖像存儲(chǔ)在全局存儲(chǔ)器中，內(nèi)核函數(shù)內(nèi)部使用寄存器變量臨時(shí)存儲(chǔ)計(jì)算 用的變量，計(jì)算過程中不訪問全局存儲(chǔ)器，每個(gè)內(nèi)核函數(shù)完成X個(gè)均勻分布視角的反投影 計(jì)算，計(jì)算完X個(gè)視角后再給全局存儲(chǔ)器的相應(yīng)像素位置賦值，以便節(jié)省反投影重建的訪 存時(shí)間，根據(jù)具體應(yīng)用環(huán)境不同，該X=N/M。N為掃描一圈采集的全部視角；M為用來(lái)并行 處理視角的線程數(shù)量，且能夠被N整除，M等于1時(shí)，對(duì)重建的像素點(diǎn)不存在"寫沖突"，可以 不使用原子加操作進(jìn)行寫入全局存儲(chǔ)器，當(dāng)M大于1時(shí)，