computed tomography ( ct ) has become one of the most popular diagnostic modalities since its invention thirty years ago . compared with 2d parallel - beam and fan - beam ct , 3d cone - beam ct system is able to achieve higher special resolution and better utilization of photons . with the rapid development of detector technology , the single detector unit is getting smaller and smaller while the number of detector units is becoming larger and larger . this means that there will be larger amount of projection data needed to be processed in 3d cone - beam ct system . paxscan2520 , a flat panel detector made by varian , has 1920  1536   detector units . the size of reconstruction image is about 512 mb for a 512 image array and is 4 gb for a 1024 image array . currently , image reconstruction speed is still a bottleneck for the development of 3d cone - beam ct . the study of fast and efficient reconstruction algorithms for large volume image and their implementation on hardware or software will have important significance both theoretically and practically [ 2 , 3 ] .                                         the graphics processing unit ( gpu ) can process volume data in parallel when working in single instruction multiple data ( simd ) mode . because of the increasing demand of computer game market and engineering design , the development of gpu has been much faster than cpu . nowadays the increasing programmability of gpu has made it possible that certain general purpose computing based on cpu can be implemented on gpu with a much faster computation speed , and general purpose gpu computing has become another hot research topic , which includes its application on ct image reconstruction [ 3 , 5 ] .                           back to 1990s , only high - end workstations , such as the sgi octane or onyx , had the level of graphics hardware necessary for ct image reconstruction . were the first to employ this hardware for the acceleration of ct reconstruction . with the fast development of the low - cost pc - based graphics hardware of similar capabilities than that of the sgi , many researchers recently have tried to implement both iterative and analytic ct image reconstruction using gpu - acceleration [ 5 , 710 ] . usually the size of graphics card memory is much less than system memory , which imposes a big constraint on the gpu - based image reconstruction . the size of the volume image array for most gpu - accelerated ct reconstruction is generally limited to 512 . they divided the target volume using a method similar to octree and proposed to reconstruction those small bricks one by one for large volume data . also presented a memory management strategy that decreases the bus transfer between main memory and gpu memory for reconstructing large volume . in this paper , we study how to partition the data to fit them into the graphics card memory . a new method of the projection data partition for large volume data is proposed . according to rotational symmetry and vertical symmetry in circular cone - beam projection layout , the method divides both projection data and reconstructed image volume into subsets . by packing the subsets of data into the rgba channels of gpu , a fast reconstruction for large data volume can be implemented . our gpu - accelerated backward - projection for fdk algorithm is introduced , and then the utilization of geometric symmetries is described . in section 3 , the partition scheme for reconstructing large data volume is presented . in section 4 , numerical experiments on various datasets are presented to evaluate the speedups with our method . in section 5 , the filtered back - projection algorithm proposed by feldkamp , davis and kress ( fdk ) for 3d volume reconstruction from circular cone - beam projections still remains one of the most widely used approach . in this algorithm , the most time - consuming part is the back - projection procedure , which has a complexity of o(n ) in the spatial domain and constitutes the bottleneck for all software solutions . so here we purposely concentrate on using gpu to accelerate the backward - projection of the reconstruction program . because the backward - projection is very similar for different ct image reconstruction algorithms , it will be easy to adapt our scheme into different reconstruction algorithms . current gpus can be used either as a graphical pipeline or as a multiprocessor chip thanks to the cuda interface from nvidia . for both options , xu and mueller have observed that an implementation of the cone beam back - projection using the graphics pipeline is 3 times faster than the one made with cuda interface . hence we use the graphics pipeline to accelerate ct reconstruction in this paper . in order to harness gpu to provide acceleration of 3d volume image reconstruction , we represent reconstructed volume as an axis - aligned stack of 2d - textured slices . the volume may be represented by three kinds of proxy geometries as shown in figure 1 . if the stack of 2d - textured slices aligned along the rotation axis yv ( figure 1(c ) ) is adopted , only one data set is enough for circular cone - beam reconstruction . otherwise , two copies of the data set should be used simultaneously in gpu memory for decreasing the inconsistent sampling rate of volume . this can cause bottlenecks when the memory bandwidth is less than the compute bandwidth , and also needs to merge the two textured slices stacks in each backward - projection loop . hence , we choose the model shown in figure 1(c ) as reconstructed volume model . figure 2 shows the geometry between the x - ray source , the reconstructed volume and the detector of cone - beam ct . the source - to - rotation center distance is d , the source - to - detector distance is d , and the rotation axis of reconstructed volume is yv . the target volume is represented as a stack of 2d slices ( textures ) perpendicularly aligned along yv axis . the key for the gpu - accelerated backward - projection is to calculate the projection positions of the vertices of every slice on the x - ray detector of cone - beam ct system . if the projection positions of the four vertices of a slice on the x - ray detector are produced , we can generate the projection coordinate of each voxel of the slice by interpolating the coordinates of the vertices in gpu rasterizer , and achieve the backward - projection from a projection image to the slice in gpu fragment shader . here is our detailed algorithm for backward - projecting a projection image to a volume slice based on gpu . as shown in figure 2 , the slice under reconstruction has four vertices v1 , v2 , v3 , and v4 , whose projection positions in detector space are p1 , p2 , p3 , and p4 , respectively . according to projective - texture mapping theory [ 14 , 15 ] , we decompose and express the full coordinate transformation from volume space to detector space as a series of matrices , as shown in formula ( 1 ) . as compared with the method presented in , formula ( 1 ) focuses on calculating the projection coordinates of any one vertex of a slice in circular cone - beam projection layout . the coordinates in detector space of the four projection positions may be computed in parallel in gpu according to this formula :  
 ( 1)eptrv = vp,[nw200nw20nh20nh200110001][2dw00002dh0000(f+n)fn2fnfn0010][10000100001d0001][cos  0sin00100sin0cos  00001][xvyvzv1]=[xpypzpwp ] . 
  in formula ( 1 ) , the coordinate of a slice vertex in volume coordinate system is expressed as a 4d homogenous vector v = ( xv , yv , zv , 1 ) . a 4  4   rotation matrix r rotates the volume coordinate system by      degrees in counter - clockwise direction . another 4  4 translation matrix t   translates the volume coordinate system a distance of d along the negative axis . the two matrices   r   and t   jointly map the vertex coordinate v from volume coordinate system ( xv  yv  zv ) into source coordinate system ( xs  ys  zs ) . a 4  4 perspective projection matrix p , determined by the source location and the detector dimensions   w    and h ,   defines a frustum for cone - beam projection . n   and   f of the matrix p denote the distances from x - ray source to the near and far clipping planes of the frustum , respectively . the matrix p implements the subsequent perspective projection that maps the frustum into a cube clip space , whose cartesian coordinates are between   1 and 1 . then a texture coordinate conversion matrix e , defined by the horizontal and vertical numbers of detector units , produces the homogeneous coordinates vp in detector space for the vertex v.                          by implementing the calculation of formula ( 1 ) in simd way in the vertex shader of gpu , the rectangle texture homogeneous coordinates of the four projection positions   p1 , p2 , p3 , and p4 in detector space are obtained . then the fragments corresponding to the voxels of slice are generated in orthographic viewing mode in gpu rasterizer , and the texture homogeneous coordinate of each fragment is produced by the linear interpolation of the texture homogeneous coordinates of p1 , p2 , p3 , and   p4 . to compensate for the perspective distortion effects , the texture coordinate ( xp , yp ) of each fragment is divided by its 4th component wp to derive correct coordinate in the fragment shader of gpu . at last , these texture coordinates are used to sample the projection image of this projection view , and the obtained sample values are accumulated into the corresponding voxels of the output texture representing the slice . note that the sample positions usually do not coincide with the detector units , the final values of the sample positions are produced by nearest - neighbor interpolation or bilinear interpolation . the above procedure is executed repeatedly until every volume slice is processed from every projection view , thus the entire reconstructed volume is updated . in circular cone - beam volume reconstruction , there are two types of geometric symmetries , which are referred to as the rotational symmetry [ 2 , 7 ] and vertical symmetry . that is , the pair of the x - ray source s1 and the voxel position v1 can be replicated by rotating it across 90 , 180 , and 270 intervals respectively to produce the other three pairs of ( s2 , v2 ) ,   ( s3 , v3 ) , and ( s4 , v4 ) . the backward - projection can be significantly speeded up by the utilization of rotational symmetry , since the geometry transform matrix described in formula ( 1 ) and sample positions in projection images are calculated only once for four symmetric projection views . there is still another kind of rotational symmetry , that is , two pairs of source and pixel positions are symmetric with respect to a diagonal line , which is also called complement symmetry . constrained by the inherent parallelism offered by the four color channels of gpu , we only utilize the 90-degree symmetry to accelerate back - projection by packing four rotational symmetric projections into four color channels . our gpu - based backward - projection algorithm using the 90-degree symmetry is as following . step 1arrange the projection images in the four rotational symmetric views of  ,     + 90 ,     + 180 , and  + 270   as one group , and pack them into the four color channels ( red / green / blue / alpha or rgba ) of a 2d - texture projtex , one projection image per channel . arrange the projection images in the four rotational symmetric views of  ,     + 90 ,     + 180 , and  + 270   as one group , and pack them into the four color channels ( red / green / blue / alpha or rgba ) of a 2d - texture projtex , one projection image per channel . step 2employ four textures slicetex1 , slicetex2 , slicetex3 and slicetex4 to save the backward - projected values for four slices , respectively . each of the four textures has individual   four color channels , and each channel is used to save the backward - projected values from one projection view . employ four textures slicetex1 , slicetex2 , slicetex3 and slicetex4 to save the backward - projected values for four slices , respectively . each of the four textures has individual   four color channels , and each channel is used to save the backward - projected values from one projection view . step 3 in gpu vertex shader , the computation described in formula ( 1 ) is only done once from projection view    for each slice by using the algorithm presented at section 2.1 , to produce the projection texture coordinates that are identical for the four symmetric projection images . in gpu vertex shader , the computation described in formula ( 1 ) is only done once from projection view    for each slice by using the algorithm presented at section 2.1 , to produce the projection texture coordinates that are identical for the four symmetric projection images . step 4 in gpu fragment shader , the projection images in the four symmetric views are backward - projected and accumulated to the four slice textures , respectively , according to the projection texture coordinates produced by the vertex shader and the succeeding rasterizer of gpu . the four slice textures are then rendered to gpu frame buffer in the same pass by the multiple render targets ( mrt ) technique of opengl . in gpu fragment shader , the projection images in the four symmetric views are backward - projected and accumulated to the four slice textures , respectively , according to the projection texture coordinates produced by the vertex shader and the succeeding rasterizer of gpu . the four slice textures are then rendered to gpu frame buffer in the same pass by the multiple render targets ( mrt ) technique of opengl . step 5the above procedures are repeated with ping - pong technique , until the projection images of all views are backward - projected and accumulated to the four slices . the effort of backward - projection from full 360 arc is reduced to one 90 arc by using the rotational symmetry . the above procedures are repeated with ping - pong technique , until the projection images of all views are backward - projected and accumulated to the four slices . the effort of backward - projection from full 360 arc is reduced to one 90 arc by using the rotational symmetry . pass , as shown in figure 4 , the texture coordinates in the g , b and a channels of the four slice textures are rotated by 90,180 , and 270 , respectively , and the data in the rgba channels of each slice texture are respectively accumulated and packed into an output texture with four channels , one slice per channel , which is then downloaded to system memory . now the four slices have been updated by the projection data from all the views on circular trajectory . the method increases the speed of downloading data by taking advantage of the 4-channel rgba parallelism , and avoids the calculation of slices accumulating in cpu . pass , as shown in figure 4 , the texture coordinates in the g , b and a channels of the four slice textures are rotated by 90,180 , and 270 , respectively , and the data in the rgba channels of each slice texture are respectively accumulated and packed into an output texture with four channels , one slice per channel , which is then downloaded to system memory . now the four slices have been updated by the projection data from all the views on circular trajectory . the method increases the speed of downloading data by taking advantage of the 4-channel rgba parallelism , and avoids the calculation of slices accumulating in cpu . step 7the above processes are repeated for every volume slice from every projection view , and then the entire reconstructed volume is updated . the above processes are repeated for every volume slice from every projection view , and then the entire reconstructed volume is updated . another type of symmetry is known as vertical symmetry in circular cone - beam projection layout . as shown in figure 5 , the vertices v1 and v2 of the reconstructed volume are vertical symmetric with respect to the central scanning plane ( zv ,   xv ) , that means when the coordinate of v1 is ( xv , yv , zv ) , the coordinate of v2 is ( xv ,   yv , zv ) . according to formula ( 1 ) , if the projection coordinate of v1 in detector space is ( xd ,   yd ) , then that of v2 is ( xd ,    yd ) in the circular cone - beam scanning case . that is , their xd coordinates stay the same , while their yd coordinates are opposite . we use the property of vertical symmetry to decrease the amount of backward - projection positions calculation . when loading a projection image into gpu memory , we read the data in the upper half of the projection image along yd axis , but read the data in the lower half in the opposite order of yd   axis , that is , fold the projection image . then we pack the two halves of the projection image into two color channels of a 2d - texture , respectively . thus the projection positions of the vertices of two vertical symmetrical slices in the projection image are identical , consequently only half of projection positions are needed to calculate for backward - projecting a projection image to the reconstructed volume . as for gpu - accelerated algorithms , the projection data should be firstly loaded into graphic card memory so as to be called by gpu , which required expensive data transfers between graphic card memory and system memory due to bandwidth limit . since the reconstruction of each slice needs the projection images from all projection views , we try to load all the projection images into graphics card memory at one time for saving data transfer time . if the amount of projection data exceeds the graphic card memory capacity , the projection data have to be partitioned into blocks to fit into the graphic card . a new partitioning scheme is employed in our program . as shown in figure 6 , the reconstructed volume is divided into several chunks , each of which is a stack of the slices of volume . the projection data for reconstructing a chunk do not require complete full sized projection images , but only the blocks contained in a rectangular shape . the height of the rectangle is greatest when the diagonal of volume slices is perpendicular to the detector plane . considering these properties , we divide projection images into the same number of blocks as volume chunks . the size of each block may be calculated by the following formulas ( 2)(4 ) . let   d   be source - to - rotation center distance , d   be source - to - detector distance , and hd    be detector height , then we can get the maximum height of the reconstructed volume : 
 ( 2)hv = d    hdd+(2/2)hd . 
  if the reconstruction volume is divided into n chunks , then the height of each chunk is hv / n , and the top coordinate and bottom coordinate of the nth chunk along the yv axis can be calculated according to the formula ( 3 ) :  
 ( 3)vtn = hv(12nn)vbn = hv(12n+1n )  n=0,1,  ,n1 . the   corresponding top projection position tn and the bottom projection position bn of the nth chunk along yd axis in detector may then be calculated by formula ( 4 ) , as shown in figure 6 : 
 ( 4)tn = d    vtnd(2/2)hvbn = d    vbnd+(2/2)hv  n=0,1,  ,n1 . 
  according to the obtained parameters tn and only one related projection data block is uploaded into graphics card memory each time for reconstructing one chunk , and backward - projected to all slices of the chunk . the actual number of partitions in the program will depend on the size of projection data and the size of graphic memory . different from the partitioning scheme simply introduced in , our method decomposes projections into blocks by utilizing vertical symmetry and rotational symmetry , and reconstructs the volume slices at the vertical symmetrical positions in the meanwhile . the partitioning scheme can avoid repeated data transfer and speed up volume reconstruction .                  in conjunction with the geometric symmetries presented in sections 2.2   and 2.3 , the cone - beam ct image reconstruction method for large data volume is summarized as follows . step 1the projection at each view is divided into blocks according to our partitioning scheme , and each block is further decomposed into two vertical symmetric subblocks for utilizing vertical symmetry , that is , an upper sub - block and a lower sub - block . the projection at each view is divided into blocks according to our partitioning scheme , and each block is further decomposed into two vertical symmetric subblocks for utilizing vertical symmetry , that is , an upper sub - block and a lower sub - block . step 2the data in the upper sub - block of current block are packed into a texture with rgba color channels every four rotational symmetric views , and the data in the lower sub - block are also packed into another texture with rgba color channels every four rotational symmetric views . all data in the current block are transferred into graphic card memory from system memory by these textures for subsequent image reconstruction . the data in the upper sub - block of current block are packed into a texture with rgba color channels every four rotational symmetric views , and the data in the lower sub - block are also packed into another texture with rgba color channels every four rotational symmetric views . all data in the current block are transferred into graphic card memory from system memory by these textures for subsequent image reconstruction . step 3according to the algorithms presented in sections 2.2   and 2.3 , four vertical symmetric slices are reconstructed in parallel each time using the current projection data block from four rotational symmetric views . the process is repeated until all slices in current data chunk are reconstructed . according to the algorithms presented in sections 2.2   and 2.3 , four vertical symmetric slices are reconstructed in parallel each time using the current projection data block from four rotational symmetric views . step 4once the image reconstruction for every four slices is achieved , they are packed into an output texture with four channels , one slice per channel , and downloaded into system memory . once the image reconstruction for every four slices is achieved , they are packed into an output texture with four channels , one slice per channel , and downloaded into system memory . step 5the above 2nd to 4th steps are executed repeatedly until every volume slice in every chunk is processed from all projection views , then the image reconstruction for entire volume is achieved . the above 2nd to 4th steps are executed repeatedly until every volume slice in every chunk is processed from all projection views , then the image reconstruction for entire volume is achieved . to test the gain of our gpu - based acceleration scheme , we have used the fdk algorithm that applies the gpu - based backward - projection to reconstruct images from computer simulated data and real mouse data acquired with a microcone - beam ct system . the pc used has a 1.83 ghz intel xeon 5120 dual - core cpu with 8 gb of system memory . the graphics card is nvidia qurdro fx4600 model with 768 mb of memory . for the simulated data , the source - to - rotation center ( sod ) is set to 1660.0 mm , the source - to - detector distance ( sdd ) is 1900.0 mm , and the size of each detector bin is 0.127 mm      0.127 mm . the micro - ct source - to - detector distance is 570.0 mm and source - to - rotation center is 390.0 mm . we have performed reconstructions for shepp - logan phantom volumes with 512 and 1024 voxels by use of the fdk algorithm with a gpu - based backward - projection . in this reconstruction , the detector array sizes are 512 and 1024 , and the numbers of projection views are 360 and 720 , respectively . the programmable pipeline of fx4600 gpu supports 32-bit float precision calculation , and our gpu - based reconstructions show the equivalent image quality as our cpu - based implementations , as shown in figure 7 . since our graphics card memory is 768 mb , the projections for reconstructing the volume with 512 voxels can be uploaded into graphics card memory at one time , and the backward - projection takes 7.27.7 seconds . while the projections in 32-bit float precision for reconstructing the volume with 1024 voxels are too large to be transferred to graphics card memory at one time . the projections need to be partitioned to fit into the graphics card memory on the basis of our partitioning scheme presented in section 3 . we decompose the reconstruction volume into 4 chunks , and correspondingly the projection data are also divided into 4 blocks . since the height of each detector unit is 0.127 mm , the height of the detector hd    is 1024  0.127 mm = 130.038 mm . according to formulas ( 2)(4 ) , we can know the maximal size of reconstruction volume hv = 108.368 mm , and the top projection coordinates and bottom projection coordinates of these 4 chunks are : t0 = 65.024 mm , b0 = 29.643 mm ; t1 = 32.512 mm , b1 = 0 ; t2 = 0 , b2 =   32.512 mm ; t3 =   29.643 mm , b3 =   65.024 mm , respectively . that is , about 279 rows of projection data at each view are needed for reconstructing the 0th chunk of volume , 256 rows for the 1th chunk , also 256 rows for the 2th chunk , and 279 rows for the 3th chunk . altogether ( 279 + 256 + 256 + 279)/1024  1.045 times projection data are needed to transfer from system memory to graphics card memory for reconstructing the whole volume . hence , as compared to the methods presented in papers [ 5 , 9 ] , the amount of transferred data is greatly reduced by our division scheme , and the backward - projection time is only 101.9104.8 seconds . we have also applied the gpu - accelerated fdk algorithm to the real mouse data acquired with a microcone - beam ct scanner . the projection size for each projection view is 1600  980 , and data were collected at a total number of 360 views . again , this system is too large for one shot reconstruction , and the projection data needs to be partitioned . the backward - projection time is about 14.515.2 seconds given by 3 partitions of the projection data . figure 8 shows the sagittal slice and the middle transverse slice of the reconstruction image . as compared to our cpu - based implementation of the fdk algorithm on the same computer , the backward - projection time in the work , we have investigated and implemented a gpu - based 3d cone - beam ct image reconstruction algorithm for large data volume , and evaluated the gpu - based implementations by use of computer - simulation data and real micro - ct data . the gpu - based implementation using geometric symmetries has speeded up the backward - projection process by about 110120 times for a circular cone - beam scan , as compared to the cpu - based implementation on the same pc . further work is in progress to apply our algorithms to the iterative image reconstruction methods of cone - beam ct .
currently , 3d cone - beam ct image reconstruction speed is still a severe limitation for clinical application . the computational power of modern graphics processing units ( gpus ) has been harnessed to provide impressive acceleration of 3d volume image reconstruction . for extra large data volume exceeding the physical graphic memory of gpu , a straightforward compromise is to divide data volume into blocks . 
 different from the conventional octree partition method , a new partition scheme is proposed in this paper . 
 this method divides both projection data and reconstructed image volume into subsets according to geometric symmetries in circular cone - beam projection layout , and a fast reconstruction for large data volume can be implemented by packing the subsets of projection data into the rgba channels of gpu , performing the reconstruction chunk by chunk and combining the individual results in the end . 
 the method is evaluated by reconstructing 3d images from computer - simulation data and real micro - ct data . 
 our results indicate that the gpu implementation can maintain original precision and speed up the reconstruction process by 110120 times for circular cone - beam scan , as compared to traditional cpu implementation .