本發明涉及醫學領域與電子技術領域的交叉，具體涉及X射線CT技術領域中的基于CUDA架構加速CT圖像重建的方法。

背景技術：
CT(ComputedTomography)重建尤其是三維重建，計算量大、耗時高，計算復雜度與被重建體數據量、投影視圖個數的乘積成正比，比如從360個投影視圖重建512張512×512大小的圖像（即5123volume）的計算復雜度為360×5123。如何提高重建速度受到越來越多的人重視，在2011年召開的第十一屆Fully3D(The11thInternationalMeetingonFullyThree-DimensionalImageReconstruction)會議論文集里有約1/4的文章涉及到三維加速重建,在其他雜志上涉及CT重建加速的文章近年來也很多。GPU的單指令多數據流（SingleInstructionMultipleData，縮寫為SIMD）處理模式為可并行地對大規模數據進行同樣的操作。由于計算機游戲和工程設計的巨大市場驅動，GPU的發展速度大大超過了CPU的發展速度，圖形流水線的高速度和高帶寬極大地提高了圖形處理能力，近年來發展起來的可編程功能為圖形處理之外的通用計算提供了高性價比的運行平臺，使得基于GPU的通用計算成為近年來的研究熱點之一。FDK重建算法于1984年由Fedlkamp等人首先提出，對CT近似重建有著重大意義，目前廣泛應用于錐形束投影重建，而且各個角度的反投影無數據交換，具有高度的并行性，因此特別適于GPU這種單指令多數據（SIMD）的流式計算架構。最早的GPGPU（GeneralPurposeGPU，即通用計算圖形處理器）開發直接使用了圖形學API編程。這種開發方式要求編程人員將數據打包成紋理，將計算任務映射為對紋理的渲染過程，用匯編或者高級著色器語言（如GLSL、Cg、HLSL）編寫shader程序，然后通過圖形學API(Direct3D、OpenGL)執行。這種“曲線救國”的方式不僅要求熟悉需要實現的計算和并行算法，還要對圖形學硬件和編程接口有深入的了解。由于開發難度大，傳統GPGPU沒有被廣泛應用。CUDA（ComputeUnifiedDeviceArchitecture，統一計算設備架構）的GPU采用了統一處理架構，可以更加有效地利用過去分布在頂點渲染器和像素渲染器的計算資源；而且引入了片內共享存儲器，支持隨機寫入（scatter）和線程間通信。

技術實現要素：
鑒于現有技術的不足，本發明旨在于提供一種快速的CT圖像重建的異步并行處理方法，針對目前錐形束重建時數據輸入、數據加權、數據濾波及反投影串行執行的瓶頸問題，提出了基于CUDA架構中GPU異步并行處理的重建方法，從而提高重建速度。本發明的技術方案具體是：通過使用兩個或者兩個以上的流，使應用程序實現任務級的并行化，進一步地說，即GPU在執行核函數的同時，還能在主機與設備之間執行復制操作。為了實現上述目的，本發明采用的技術方案如下：一種基于CUDA架構加速CT圖像重建的方法，包括數據輸入模塊，基于GPU的CT數據加權濾波模塊，基于GPU的CT圖像重建反投影模塊，以及數據輸出模塊，其特征在于，所述方法包括以下步驟：（1）從掃描的X射線強度數據獲得投影數據，經過預處理后，由CPU讀入到內存中；（2）應用程序實現任務級的并行化，通過使用兩個或兩個以上的流，使GPU在執行核函數的同時，能在主機與設備之間執行復制操作。需要說明的是，所述CT數據加權濾波模塊在GPU中執行，為每個待加權濾波元素分配至GPU中的單獨線程來執行，其中，所述線程分配過程如下：根據GPU的特性設置每個線程塊的尺寸；根據補零后投影數據的水平長度和垂直長度設置所述線程塊的個數；按照所述線程塊設置執行內核程序。需要說明的是，所述CT圖像重建反投影模塊在GPU中執行，為每個待重建像素分配至GPU中的單獨線程來執行，重建所需濾波后的數據存儲在GPU的紋理內存中，其中線程分配過程如下：根據GPU的特性設置每個線程塊的尺寸；根據待重建圖像的尺寸設置所述線程塊的個數；按照所述線程塊設置執行內核程序。需要說明的是，所述投影數據使用所述基于GPU的CT數據加權濾波模塊，以及所述基于GPU的CT反投影模塊采用濾波反投影算法獲得重建體；其中，所述基于GPU的CT數據濾波模塊，將投影數據首先在GPU上進行加權，再通過GPU上FFT變換到頻域，頻域濾波后通過GPU上的逆FFT獲得濾波后的數據。需要說明的是，所述基于GPU的圖像重建反投影模塊在GPU上實現紋理綁定,將顯存中的數據與紋理參照系相關聯，并進行紋理拾取操作。作為一種優選的方案，紋理緩存中的數據可以被重復利用，而且一次拾取坐標對于位置附近的幾個像元，提高一定局部性的訪存效率。需要說明的是，所述兩個或兩個以上的流處理數據互不影響。需要說明的是，其特征在于，數據從內存復制到顯存、GPU上的投影數據加權操作、GPU上的投影數據濾波操作及GPU上的CT圖像重建的反投影操作為異步并行執行。需要說明的是，輸入數據存儲為無符號短整型；GPU的CT加權濾波數據、GPU的CT圖像重建反投影數據及CPU的輸出數據存儲成32位浮點格式。本發明有益效果在于，采用了異步并行的執行方法，明顯提高了CT圖像的重建速度。附圖說明圖1為平行探測器錐束掃描幾何結構圖；圖2為本發明的方法流程示意圖，其中a為開始部分流程圖；b為本循環部分流程圖；c為結束部分流程圖。具體實施方式下面將結合附圖對本發明作進一步的描述。如圖1所示，為平板探測器錐束掃描幾何結構，射線源到旋轉中心的距離為R，射線源到探測器距離為D，扇角為γ，錐角為τ，稱射線源到探測器中心且與探測器垂直的射線為中心射線，FDK算法重建公式為：其中gI(u,v,λ)代表投影數據，λ為投影角度。FDK算法實現步驟為：（1）加權濾波：（2）加權反投影：其中，U(x,y,λ)＝R+xcosλ+ysinλ需要說明的是，錐形束重建GPU模塊主要包括兩部分：CT投影數據加權濾波模塊和CT圖像的反投影模塊。假設通過X射線平板探測器獲得K個角度的二維投影數據分別為p0p1......pk-1，每一副投影由U×V個像素組成，將從中重建L×W×H的體數據F。1、基于CUDA的CT數據加權濾波模塊（1）生成加權函數并保存在顯存數組d_weight[U][V]中；（2）利用FFT變換將二維投影數據pn(0≤n≤k-1)轉換到頻域。由于對投影數據需要進行一維頻域濾波，在實現投影數據的濾波之前，需要生成窗函數，并將其做FFT變換。GPU上數據的FFT變換分以下幾步實現：第一步，依次復制內存中各個角度的二維投影數據pn(0≤n≤k-1)到顯存，記為d_inData[V][U]；第二步，對投影數據在水平方向上補零，補零后的數據長度為U′;需要說明的是，此時需要考慮以下三個因素：（a）避免干涉效應，補零的最小數目為探測器長度減1（即：U-1）；（b）實現快速FFT變換，補零后的長度應為2的整數次冪；（c）實函數經過FFT變換后在頻域上為偶函數。（3）在顯存中開辟一個二維數組d_data[V][U′],并將每張投影數據d_inData[V][U]中的每一個元素與d_weight[V][U]中的對應元素相乘，做加權操作，每行的末尾補零。intx=__mul24(blockDim.x,blockIdx.x)+threadIdx.x;inty=__mul24(blockDim.y,blockIdx.y)+threadIdx.y;if(x<U&&y<V)d_data[y][x]=d_inData[y][x]*d_weight[y][x];（4）設置CUDA濾波時參數，主要步驟為：第一步，根據CUDA的特性設置每個線程塊(Block)的尺寸；根據補零后投影數據的水平長度U′及垂直長度V設置線程塊(Block)的個數。第二步，為FFT變換和逆FFT變換分別創建FFT句柄和一維FFT句柄plan。cufftHandleplanF,planI;cufftPlan1d(&planF,U′,CUFFT_R2C,V);cufftPlan1d(&planI,U′,CUFFT_C2R,V);第三步，將FFT句柄plan與CUDA流相關聯cufftSetStream(planF,stream1);cufftSetStream(planI,stream2);（5）將投影數據進行原地（inplace）FFT變換，并與濾波窗函數的頻域值在相應位置上進行點乘，得到濾波后的數據。cufftExecR2C(planF,(cufftReal*)d_data,(cufftComplex*d_data)（6）將濾波后的數據進行原地（inplace）逆FFT變換，此時數據的水平大小仍為U′，垂直大小為V，數據仍然存儲在d_data[V][U′]中。cufftExecC2R(planI,(cufftComplex*d_data,(cufftReal*)d_data);2、基于GPU的CT數據反投影模塊在主機端聲明需要綁定到紋理的CUDA數組并設置好紋理參照系，然后將紋理參照系與CUDA數組綁定建立紋理坐標系，之后就可在內核中通過紋理拾取函數訪問紋理存儲器。具體步驟為：第一步，聲明紋理參照系，紋理參照系通過一個作用范圍為全文件的texture型變量聲明，并且必須在編譯前顯示聲明texture<float,2,cudaReadModeElementType>texRef;第二步，設置運行時紋理參考系屬性texRef1.addressMode[0]=cudaAddressModeWrap;texRef1.addressMode[1]=cudaAddressModeWrap;texRef1.filterMode=cudaFilterModeLinear;texRef1.normalized=false;第三步，根據GPU的特性設置每個線程塊(Block)的尺寸；根據待重建圖像的尺寸設置線程塊(Block)的個數；第四步，根據探測器水平大小(U)和垂直大小(V),聲明CUDA數組，并分配空間cudaChannelFormatDescchannelDesc=cudaCreateChannelDesc(32,0,0,0,cudaChannelFormatKindFloat);cudaArray*cuArray;cudaMallocArray(&cuArray,&channelDesc,U，V);第五步，將濾波后的投影數據d_data[V][U′]復制到CUDA數組cudaArray中第六步，紋理綁定，將顯存中的數據與紋理參照系相關聯的操作。cudaBindTextureToArray(texRef,cuArray,channelDesc);第七步，紋理拾取，采用紋理坐標對紋理存儲器進行訪問，即可得到體數據的值。獲取重建體各個像素點的位置：intx=__mul24(blockDim.x,blockIdx.x)+threadIdx.x;inty=__mul24(blockDim.y,blockIdx.y)+threadIdx.y;intz=__mul24(blockDim.z,blockIdx.z)+threadIdx.z;計算各像素點映射到探測器的位置,在此假設水平方向為point_h,垂直方向為point_v。則某個角度的投影數據p在該像素點的貢獻值為：tex2D(texRef,point_h+0.5,point_v+0.5)；下面結合具體實施例進一步的描述本發明。如圖2所示，利用本發明的基于異步并行處理的FDK重建方法進行重建，GPU使用NVIDIA的GeForceGT640。需要說明的是，投影數據大小為512*512*480，待重建體的大小為512*512*512。1、從外部設備（如硬盤）讀取投影數據到內存，由于投影數據存儲空間遠小于內存，待重建體存儲空間小于顯存的一半，因此投影數據全部讀取到內存，并直接生成重建體。如果投影數據過大，可以考慮分塊輸入；如果重建體過大，可以考慮分塊重建。此處假設內存中的投影數據為h_indata。2、初始化四個流對象cudaStream_tstream[4];for(inti=0;i<4;i++)cudaStreamCreate(&steam[i]);3、四個流進行異步并行操作，提高重建速度，其中：（1）第0個流從h_indata取出一張投影數據，并將其復制到顯存；（2）在第0個流做加權操作的同時，第1個流從h_indata取出下一張投影數據，并將其復制到顯存；（3）在第0個流做濾波操作的同時，第1個流進行加權操作，第2個流從h_indata取出下一張投影數據，并將其復制到顯存；（4）在第0個流做反投影操作的同時，在第1個流做濾波操作的同時，第2個流進行加權操作，第3個流從h_indata取出下一張投影數據，并將其復制到顯存；重復執行（1）～（4）操作，直到所有投影數據全部被讀入。通過使用stream，可以將數據讀取時間、加權時間和濾波時間部分隱藏，從而提高程序執行效率。對于本領域的技術人員來說，可根據以上描述的技術方案以及構思，做出其它各種相應的改變以及變形，而所有的這些改變以及變形都應該屬于本發明權利要求的保護范圍之內。