CN103077547A - CT (computerized tomography) on-line reconstruction and real-time visualization method based on CUDA (compute unified device architecture) - Google Patents

Abstract

The invention provides a CT (computerized tomography) on-line reconstruction and visualization method based on a CUDA (compute unified device architecture). The method comprises the following steps of obtaining projection data and pre-treating the obtained projection data; realizing FDK (Feldkamp, Davis and Kress) weighted filter process by a CPU (central processing unit); accelerating CUDA to realize FDK weighted back protection; and accelerating CUDA to realize volume rendering. The method provided by the invention can be used for implementing on-line reconstruction for cone beam CT data, so that a cone beam CT system can reconstruct a CT three-dimensional image while sampling the projection data, and an integrated CT image can be obtained finally along with the completion of the projection data, thereby realizing real-time on-line feedback.

Description

Rebuild and real time visualized method based on the CT of CUDA framework is online

Technical field

The present invention relates to the online method of rebuilding of X ray CT image, particularly relate to the method for the online reconstruction of a kind of pyramidal CT image based on the CUDA framework and real-time visual.

Background technology

Along with the development of hardware, CT scan speed is improved, and the image reconstruction time becomes the bottleneck of application greater than sweep time in actual applications.CT rebuilds large, the consuming time height of calculated amount, and computation complexity is directly proportional with the product of rebuilt volume data amount, projection view number, how to improve reconstruction speed and is subject to increasing people's attention.For the three-dimensional reconstruction acceleration problem, more international scientific research institutions and company have done fruitful research and development, acceleration can set about reducing from algorithm itself time complexity of algorithm, also can develop the algorithm that is fit to some specific hardware, such as rebuilding based on CPU, rebuilding, rebuild, rebuild etc. based on Cell based on FPGA based on GPU.

FDK rebuilds:

Classical FDK algorithm is by Feldkamp, and Davis, and Kress1984 proposition, FDK algorithm have done one and have been similar to, thereby obtained the three-dimensional reconstruction algorithm from two-dimentional classic algorithm.Because it has kept filtered back projection's structure, efficient stable, FDK algorithm and improvement algorithm thereof are the main flows of using all the time.The FDK efficiency of algorithm is high, can obtain in the less situation of cone angle and rebuild preferably effect.Also have at present several improved FDK algorithms to improve the cone angle problem, the FDK algorithm mainly contains two large steps:

The first step, weighted filtering:

(1) ${\tilde{p}}^C(\mathrm{\θ},a,b)=\left(\frac{R^2}{\sqrt{R^2+a^2+b^2}}{p}^C(\mathrm{\θ},a,b)\right)*g^P\left(a\right)---\left(1\right)$

P wherein ^C(θ, a, b) is data for projection,

The combined influence of fan angle and cone angle, g ^P(a) be ramp function.

Second step, weighted back projection:

$f_{\mathrm{FDK}}(x,y,z)=\underset{\mathrm{\θ}}{\mathrm{\Σ}}\frac{R^2}{U(x,y,\mathrm{\θ})}\mathrm{Interplation}\left({\tilde{p}}^C(\mathrm{\θ},a(x,y,\mathrm{\θ}),b(x,y,z,\mathrm{\θ}))\right)---\left(2\right)$

Wherein

U(x，y，θ)＝R+xcos(θ)+ysin(θ) (3)

$a(x,y,\mathrm{\θ})=R\frac{-x\sin \left(\mathrm{\θ}\right)+y\cos \left(\mathrm{\θ}\right)}{R+x\cos \left(\mathrm{\θ}\right)+y\sin \left(\mathrm{\θ}\right)}---\left(4\right)$

$b(x,y,z,\mathrm{\θ})=z\frac{R}{R+x\cos \left(\mathrm{\θ}\right)+y\sin \left(\mathrm{\θ}\right)}---\left(5\right)$

Because most applications a (x, y, θ) andb (x, y, z, θ) is non-integer, usually need interpolation.Because most applications a (c, y, θ) andb (x, y, z, θ) is non-integer, usually need interpolation.

Visual:

Volume rendering algorithm is of paramount importance a kind of in the visualized algorithm, and is well-known can produce high-quality and drawing result true to nature.The research of volume drawing is come across eighties of last century end of the eighties, and traditional volume rendering algorithm can be divided into three large classes: the rendering algorithm of image space, the rendering algorithm of object space and image and object space blend rendering algorithm.Wherein, common classical volume rendering algorithm has volume data ray cast method, Splatting algorithm and Shear Warp algorithm etc.The ray cast method is launched a light from each pixel of projection plane, passes the three-dimensional data field, and carries out the mixing of sampled point light attribute according to the order of front-to-back, finally obtains two-dimensional projection image.For certain the sampled point S in the three-dimensional data field, can adopt the method interpolation of arest neighbors, three linearities or splines to go out its light attribute.

Graphics pipeline based on GPU accelerates:

Mapping between flat panel detector and the section to be rebuild is perspective transform, this conversion is quite common in computer graphics, Brian Cabral, Deng the people at paper Accelerated volume renderingandtomographic reconstruction using texture mapping hardware.InProceedingsof the 1994symposium on Volume visualization, pages 91-98, Tysons Corner, Virginia, United States, 1994.ACM. at first use projective textures mapping (projectivetexture mapping) to finish CT at graphics workstation to rebuild and accelerate to accelerate with volume drawing, because the similarity of the back projection in the CT process of reconstruction and projective textures mapping, the most of article that uses subsequently GPU to accelerate the CT reconstruction has all used the projective textures mapping, Xu F.Xu and K.Mueller is at paper Real-time 3d computed tomographic reconstruc-tion usingcommodity graphics hardware.Physics in Medicine andBiology, 52 (12): 3405-3419,2007. middle use GeForce8800GTX video card accelerates the CT reconstruction and has obtained better effects, the GPU universal computer model of these method utilizations is actually the graphics process api function based on the video card bottom, such as OpenGL etc., this thinking itself is sought something far and wide when it is within reach, realize comparatively complicated, this accelerated method versatility is poor simultaneously, and is difficult to grasp concerning the personage of non-graphics specialty.

The CUDA framework accelerates:

In June, 2007, NVIDIA has released CUDA (compute unified devicearchitecture, unified calculation equipment framework).CUDA is a kind of with the software and hardware architecture of GPU as data parallel equipment.From then on, a large amount of parallel computations has had direct and succinct solution finally.Up to the present, the version of CUDA has carried out four versions and has promoted, and function is constantly strengthened and be perfect.

CUDA need to be by graphics API, and has adopted than the class C language that is easier to grasp and developed.The developer can be transitioned into GPU from CPU more stably from the C language of being familiar with, and needn't relearn grammer.Certainly, develop high performance GPU general-purpose computations program, the developer still needs to be grasped the knowledge of parallel algorithm and GPU framework aspect.

Compare with GPU in the past, support the GPU of CUDA at framework significantly improvement to be arranged.These two improvement make the CUDA framework be more applicable for the GPU general-purpose computations.The one, adopted unified processing framework, can more effectively utilize the computational resource that is distributed in summit renderer and pixel rendering device over; The 2nd, introduced shared storage in the sheet, support random writing and inter-thread communication.At present the computing power of GPU has far surpassed CPU, and video card capabilities still promotes with three times Moore's Law, and performance was doubled in namely per 6 months.

The CUDA programming model: as main frame, GPU can exist a main frame and several equipment as coprocessor or equipment to the CUDA programming model in a system with CPU.In this model, CPU and GPU collaborative work, Each performs its own functions.CPU is responsible for carrying out the strong issued transaction of logicality and serial computing, and GPU then is absorbed in the parallel processing task of carrying out the height threading.CPU, GPU have separate memory address space separately: the internal memory of host side and the video memory of equipment end.

The software stack of CUDA software architecture: CUDA consists of by following three layers: CUDA Library, CUDA runtime API, CUDA driver API.The core of CUDA is CUDA C language, and it comprises expands collection and a run-time library to the minimum of C language.In a program, can only use CUDA when operation API or CUDA to drive a kind of among the API, can not mix use.

The CUDA memory model: each thread has privately owned memory register and the local storage of oneself among the CUDA, and each thread block has a shared storage.At last, all threads can be accessed the same global storage among the grid.Also have two kinds of ROM (read-only memory): constant storage and texture storage device.

Summary of the invention

The purpose of this invention is to provide a kind of CUDA of utilization framework and realize the method for the online reconstruction of CT and real-time visual, so that the user can observe present reconstructed results in the CT data acquisition, thereby in time operate.

A kind of based on online the reconstruction and visualization method of unified calculation equipment CUDA framework CT, comprising:

Obtain data for projection and the data for projection that obtains is carried out pre-service;

CPU realizes FDK weighted filtering process;

CUDA accelerates to realize the FDK weighted back projection;

CUDA accelerates to realize volume drawing.

The online reconstruction and real time visualized method of CT based on the CUDA framework of the present invention, can implement online the reconstruction to the Cone-Beam CT data, so that cone-beam CT system reconstructs the CT 3-D view in recording projection data, along with improving of data for projection finally obtains complete CT image, thereby realize real-time online feedback.

Description of drawings

Fig. 1 is the software architecture of CUDA framework;

Fig. 2 is the memory model of CUDA framework;

Fig. 3 is overall operation process flow diagram of the present invention;

Fig. 4 is obtaining and pretreated operating process of data for projection;

Fig. 5 is mouse claw experimental result;

Fig. 6 is mouse whole body experimental result.

Embodiment

Basic thought of the present invention is by CPU-GPU Heterogeneous Computing pattern, so that CPU and GPU collaborative work when gathering the original projection data, are rebuild online also real-time visual according to the data that collected, thereby realized the Real-time Feedback that CT rebuilds.

Describe method of the present invention in detail below in conjunction with accompanying drawing.Its process flow diagram of a kind of specific implementation of the present invention mainly comprises three steps as shown in Figure 3: obtain and pre-service, the CPU of data for projection realizes that FDK weighted filtering process, GPU realize that FDK weighted back projection process and GPU realize the volume drawing process.The below is introduced one by one.

Step 1: the obtaining and pre-service of data for projection

The acquisition quality of data for projection and speed dependent be in employed detector and capture card, and we carry out three-dimensional reconstruction after by preprocessing process, rebuild pseudo-shadow thereby reduce.Main data pre-service content comprises: the details in a play not acted out on stage, but told through dialogues removal of detector, explorer response Concordance, and the detector compensating bad point etc.Detector is through x-ray bombardment time t, and its output valve will be comprised of three parts: direct current biasing amount 1O, detector dark current output quantity Id and the output quantity Is that is caused by X-ray of detector.Ideally, Id and Is can increase with t linearity integral time, and the slope that changes under identical experiment condition is fixed.

As electron device, even X-ray detector is not having to have certain output in the x-ray bombardment situation yet, i.e. the direct current biasing IO of detector, and detector dark current Id causes that by the electronic noise of device the temperature of detector is lower, dark current is less.This claim direct current biasing and dark current and for the details in a play not acted out on stage, but told through dialogues data.In order to remove the impact of details in a play not acted out on stage, but told through dialogues data, before carrying out formal data for projection collection, gather first multiframe details in a play not acted out on stage, but told through dialogues data, then to average to reduce noise fluctuations and cause impact, the scanning projection data of subsequent acquisition deduct these details in a play not acted out on stage, but told through dialogues data just can realize details in a play not acted out on stage, but told through dialogues.

Because X ray at the heterogeneity of detector surface distribution and the non_uniform response between each pixel of detector, can cause the output valve of each pixel of detector under identical illuminate condition different, we are referred to as the explorer response inconsistency., need to before formal recording projection data, gather first one group of flat field data, and these flat field data are done the details in a play not acted out on stage, but told through dialogues Transformatin for this reason, scan-data is carried out flat field correction.Because the manufacture craft problem of integrated circuit, the CMOS chip has the pixel of some response abnormalities, is called bad point.Before carrying out the CT reconstruction, need to compensate these bad points.In the present invention, dead pixel compensation method is the value of using the picture element interpolation replacement bad point of surrounding normal.During without pre-service, the inconsistency of detector is obvious, and has bad point to exist, and after the data pre-service, it is very even that image becomes, and bad point has also obtained compensation.Process flow diagram sees Fig. 4 for details.

Step 2:CPU realizes FDK weighted filtering process

Can find out that from overall calculation flow process of the present invention the calculated amount of every width of cloth data for projection is very large, in order to satisfy online and real-time requirement, as far as possible so that CPU and the work of GPU highly-parallel.Among the present invention, we use CPU to finish weighted filtering process in the FDK algorithm,

(1) weight phase carries out the pointwise weighting to data for projection.In order to finish as early as possible weighting, can in advance weighted value be stored in the array, directly call during calculating.In the time of the loss dirigibility, saved computing time

(2) filtering stage, not direct convolution among the present invention, but by data for projection being carried out the zero padding operation, carry out line by line Fourier transform, then every delegation of data for projection and the filter function behind the Fourier transform carry out pointwise and multiply each other, and at last data for projection are carried out the Fu Shi inverse transformation line by line.The present invention uses the filtering in CPU of FFTW storehouse.

CPU will with process after data for projection import the GPU global storage into after, do not wait for that GPU calculates to finish, directly return and carry out next width of cloth data acquisition and process, CPU and GPU by the global storage variable realize signal synchronously, avoid occurring conflict.

Step 3:GPU realizes FDK weighted back projection process

The object space of rebuilding is successively rebuild, and the present invention uses the CUDA framework to call this weighted back projection process of concurrent execution of Kernel function, and Thread Count is determined by the gridDim parameter in its Kernel function and blockDim parameter.The present invention is fixedly installed blockDim and is 16*16, therefore has simultaneously 256 threads to carry out simultaneously among a block, and a SM comprises three block at least simultaneously, has effectively reduced delay.Simultaneously will use continually trigonometric function in back projection's process, it is very consuming time to use GPU to calculate trigonometric function, therefore can calculate first the value of trigonometric function before back projection, deposits in the geometric parameter array, directly calls according to the geometric position during calculating to get final product.This geometric parameter array can leave in the constant storage, and constant storage is read-only address space, and it is positioned at video memory, accelerates but have buffer memory.When from the same data in the thread accesses constant storage of same half-warp, if cache hit occurs, so only need one-period just can obtain data.Each SM has the constant storage buffer memory of 8KB, and it is read-only, does not have the buffer consistency problem.CT for the magnanimity data for projection rebuilds, when rebuilding i section, needed data for projection required storage space in the process of back projection may surpass the greatest physical capacity of video memory, therefore the present invention's each section of rebuilding successively object, according to the required data for projection layer of each section of trying to achieve, rebuild in batches and write back disk.

Use the CUDA framework, in global storage, open up the data after enough space storage CPU process, realize FDK rebuild in the calculating of weighted back projection, result of calculation is added in the pre-assigned CUDA array.In this process, CPU can not wait for that GPU returns, but directly begins to gather next width of cloth data and carry out pre-service.

Step 4:GPU realizes the volume drawing process

Utilize CUDA to realize light projecting algorithm, the calculating of every light is equivalent to call a thread.The present invention carries out sampling on volume data traversal, the radiation direction and the accumulation calculating of color value by stating a kernel function.Volume data is arranged in global storage after step 3, at first be set in gridDim parameter and the blockDim parameter that needs call in the Kernel function according to picture size.The present invention is made as the size of image texture data and is that 512*512*512, blockDim are 16*16, namely comprises 256 concurrent threads in each block piece.Next calculates bounding box and reduces unnecessary calculating, sets first a bounding box, and the light that only passes bounding box just might pass through object, just is necessary to calculate.Then determine coordinate and the direction of light, by the traversal of circulation realization to the volume data on the light, sampling, colouration, illumination calculation is also accumulated calculating.At last, realize the mixing of current light up-sampling point color, obtain the color value of each pixel on the image.In the process of mixing, by drawing the texture polygon of all generations from backward front order, obtain the drawing result of sub-block.Mix (Blend) computing in the following way:.

C _dst＝(1-α _src)C _dst+α _srcC _src

Wherein, Csrc is current drafting color, and Cdst is the color in the frame buffer, and α src is the alpha component (being equivalent to opacity) of current drafting color.

Operation result

We have realized above-mentioned algorithm with CUDA C, with checking algorithm complexity that the present invention is carried and practicality.Fig. 5, Fig. 6 are two examples that method of the present invention is used for the practical medical image, the size of this data for projection is 500*2340*2234 (float data), rebuilding size is 512*512*512 (float data), test platform is: 2.27GHz Intel Xeon double-core CPU, the 16GB internal memory, Tesla C1060GPU; Development environment is: Visual Studio 2010, cuda 4.0runtime API.For the single width data for projection, computing time is as shown in the table, notices that CPU-GPU is in heterogeneous schemas, and the single width total time is that GPU expends time in and the maximal value of CPU in expending time in, and adds data write GPU from CPU time.

Step	Time
		Cpu data gathers and pre-service (simulation)	10ms
CPU realizes FDK weighted filtering process	18ms
		Cpu data writes GPU	/ms
GPU realizes FDK weighted back projection process	26ms
		GPU realizes the volume drawing process	34ms

Claims (12)

1. one kind based on online the reconstruction and visualization method of unified calculation equipment CUDA framework CT, comprising:

Obtain data for projection and the data for projection that obtains is carried out pre-service;

CPU realizes FDK weighted filtering process;

CUDA accelerates to realize the FDK weighted back projection;

CUDA accelerates to realize volume drawing.

2. such as right 1 described method, it is characterized in that described pre-service comprises:

Before formal recording projection data, gather one group of flat field data and multiframe details in a play not acted out on stage, but told through dialogues data.

3. method as claimed in claim 2 is characterized in that multiframe details in a play not acted out on stage, but told through dialogues data are averaged, to reduce the impact of noise fluctuations.

4. method as claimed in claim 2 is characterized in that described flat field data are done the details in a play not acted out on stage, but told through dialogues Transformatin, scan-data is carried out flat field process.

5. method as claimed in claim 4 is characterized in that the bad point in the image is compensated, and wherein, uses the value of the picture element interpolation replacement bad point of surrounding normal.

6. the method for claim 1 is characterized in that described weighting comprises:

In advance weighted value is stored in the array, data for projection is carried out the pointwise weighting.

7. the method for claim 1 is characterized in that filtering comprises:

By data for projection being carried out the zero padding operation, carry out line by line Fourier transform;

Every delegation of data for projection and the filter function behind the Fourier transform carry out pointwise and multiply each other;

Data for projection is carried out the Fu Shi inverse transformation line by line.

8. the method for claim 1 is characterized in that using the Kernel function to carry out described weighted back projection.

9. method as claimed in claim 8 is characterized in that the blockDim in the Kernel function is fixedly installed and is 16*16, and a SM comprises three block at least simultaneously.

10. the method for claim 1 is characterized in that having calculated first the value of trigonometric function before back projection, deposit in the geometric parameter array of constant storage.

11. the method for claim 1 is characterized in that described weighted back projection comprises:

CT for the magnanimity data for projection rebuilds, and rebuilds successively each section of object, and disk is rebuild in batches and write back to the data for projection layer required according to each section.

12. the method for claim 1, it is characterized in that described CPU handles the single width data for projection after, do not wait for the calculating of GPU, directly return the collection that continues next width of cloth, described CPU and GPU realize that by global variable signal is synchronous.

Images