however , due to the high computational cost required for 
 reconstruction , its real - time imaging applications   remain 
 challenging .
bolus - chasing computed tomography ( ct ) angiography is 
 a primary example which demands real - time ct feedback . to address this problem ,
for example , xtrillion ( by terarecon , 
 inc . ) uses an application - specific pci card , while mercury 
 computer systems relies on blade - based linux clusters . however , 
 the specialized hardware is expensive and unsuitable for general 
 purpose applications .
alternatively , efforts are made using 
 graphic cards [ 2 , 3 ] , since the main operation for commercial 
 ct reconstruction is backprojection , similar to texture mapping in 
 computer graphics .
although graphics cards are highly 
 optimized , they do not support floating - point calculations .
that 
 the latest graphics cards can implement virtual floating - point 
 calculations [ 3 , 5 ] , they do not support full 32 bits 
 floating calculations .
another bottle - neck is that the graphic 
 cards require data exchange between cpu and gpu . in this paper ,
a multi - core pc - based acceleration scheme is 
 proposed for filtered - backprojection-(fbp- ) based image 
 reconstruction .
third , the single - instruction multiple - data 
 ( simd ) technique is employed for data - level parallel processing . 

fourth , the multithreading programing is done to take advantage of 
 multi - core processors , realizing the true parallel computation 
 capability . finally , an intel c++ complier is used to optimize the 
 code for intel processors .
the ct reconstruction algorithm is 
 overviewed , and then each of our acceleration techniques is 
 described . in section 3 , numerical experiments on 
 various datasets and different pcs are presented to 
 evaluate the speedups with our scheme and the conventional 
 implementation . in section 4 , relevant issues and 
 research directions are discussed .
the most popular multislice ct reconstruction methods remain data 
 rebinning - based fan - beam reconstruction filtered backprojection 
 ( fbp ) algorithms .
therefore , our work is focused on the typical 
 fan - beam fbp algorithm .
note that the application of our scheme is 
 not limited to the fan - beam case , because it can also be applied 
 to accelerate the latest approximate cone - beam algorithms 
 [ 68 ] , which can be treated as generalized fan - beam 
 reconstruction algorithms . in a typical ct setting
, the data acquisition system ( an x - ray 
 source and a detector assembly ) is rotated rapidly in the gantry 
 while the patient on a table is translated into the gantry 
 opening .
this process is illustrated in figure 1 . 
 because the multi - row detector arrays span a very small cone 
 angle , acquired helical scan data are usually rebinned into a 
 series of virtual circular scan data for reconstruction of a stack 
 of images .
here we assume a method from , in 
 which the virtual fan - beam projection data are calculated 
 according to the following formula : 
 ( 1)pc  ( ,,z0)=ph(,,z())zbza  + zb  + ph(,+2,z(+2))zaza  + zb  , 
 where  =  + k 
  2  , k  
 n , so that za 

 z0 
 zb , pc 
  denotes virtual circular scan data , ph 
  denotes acquired helical projection , 
 za and zb 
  are the distances from projections a and 
 b to the virtual circular plane , respectively , in 
 figure 2 ,  and  are the projection 
 angles shown in figure 3 . after transforming helical projection data 
 ph to circular fan - beam projection 
 data pc ,
the conventional fan - beam 
 reconstruction algorithm   can be used . as the rebinning 
 cost is insignificant ,
our optimization targets the reconstruction 
 process : 
 ( 2a)f(x , y , z0)=021l2pc , f(0,,z0)d , 
 ( 2b)pc , f(,,z0)=pc(,,z0)dcos  22sin2h( ) , 
 ( 2c)l = d2+r22drcos( ) , 
 ( 2d)0  = arcsinrsin( )  l , 
 where f is an object function to be 
 reconstructed , pc , 
 f are the filtered projection data , 
 d is the distance from the source to the center 
 of rotation , and h( ) is the ramp filter 
 .
while the inner convolution is the filtration process , 
 the outer integration is the most time - consuming backprojection 
 process , as shown in figure 4(a ) . 
 
since the backprojection is the bottleneck , let us analyze the 
 backprojection process as shown in algorithm 1 . 
 clearly , a large part of the computational cost is due to the 
 inner loop that calculates   0 , 
 1/l , interpolation coefficients , and 
 accumulates the incremental contributions to the final 
 reconstruction . in the following ,
we show how the backprojection 
 can be speeded up using various techniques . for our circular fan - beam reconstruction ,
two types of symmetries 
 are available , which are referred to as the right - angle symmetry 
 and complement symmetry .
the right - angle symmetry , or 90-degree 
 symmetry , is shown in figure 5 .
that is , a new pair of 
 source and pixel positions is obtained by applying a 90-degree 
 rotation to a current pair of source and pixel positions .
the 
 resultant 4 pairs of source and pixel positions share the same 
 1/l and 0 , which 
 can be calculated from ( 3 ) and ( 4 ) , respectively . 
 as the interpolation coefficients required by the backprojection 
 are determined by   0 , they are the same as 
 well .
therefore , for the four sets under consideration , the 
 calculations of these parameters need to be done only once : 
 ( 3)lset1  = d2+r22drcos  ( )=d2+(xp2+yp2)2drcos(tan1(xp , yp))=d2+(yp2+xp2)2drcos(+2tan1(yp,xp))=  lset2  = d2+(xp2+yp2)2drcos(+tan1(xp,yp))=  lset3=d2+(yp2+xp2)2drcos(+32tan1(yp , xp))=  lset4 , 
 ( 4)0,set1=arcsinrsin(tan1(xp , yp))l = arcsinrsin(tan1(yp , xp)/2)l=0,set2=arcsinrsin(tan1(xp , yp))l=0,set3=arcsinrsin(tan1(yp , xp)3/2)l  = 0set4 , 
 where ( xp , 
 yp ) denotes the pixel 
 in the first quadrant .
the following two requirements , which   are 
 usually satisfied in practice , are necessary to use the 
 right - angle symmetry .
the first requirement is that the projection 
 data must be available at the involved four angles .
namely , the 
 number of projections in a full scan must be divisible by 4 , which 
 is reasonable for current medical ct scanners . for instance
, a 
 somatom system generates 1160 projections per turn , while a 
 lightspeed scanner produces 984 projections per turn .
the other 
 requirement is that the reconstruction region must be symmetric 
 about the x- and y - axes , such 
 as a square or a circle in the clinical imaging situation . 
 
the second type of symmetry is the complement symmetry , as shown 
 in figure 6 . here
, a pair of source and pixel 
 positions complements the other pair of source and pixel positions 
 if they are symmetric with respect to a diagonal line ( e.g. , 
 y = x ) . for these 2 pairs of 
 source and pixel positions , l s are the 
 same , while   0 for one has the opposite 
 sign of that for the other , as shown by ( 5 ) and 
 ( 6 ) , respectively , 
 ( 5)lset1  = d2+r22drcos  ( )=d2+(xp2+yp2)2drcos(tan1(xp , yp))=d2+(yp2+xp2)2drcos(2tan1(yp , xp))=  lset2c , 
 ( 6)0,set1=arcsinrsin(tan1(xp , yp))l = arcsinrsin(tan1(yp , xp)(/2))l=0,set1c . 
 therefore , such a symmetry can also be used to reduce the 
 computational cost .
the requirements for use of the complement 
 symmetry are the same as those for the right - angle symmetry . 
 
using these two types of symmetries , the backprojection can be 
 significantly speeded up , since only one set of parameters needs 
 to be calculated for the eight sets .
the implementation of the 
 backprojection is accordingly modified , as shown in 
 algorithm 2 .
note that after the calculation of 
 0 and l once for 8 pairs of source and detector 
 positions , 8 filtered projection values are put to 8-pixel 
 positions together in the inner - loop .
to evaluate the computational complexity , the time for cpu to 
 access data must be considered , especially for the ct 
 reconstruction process because the backprojection requires 
 frequent visits to a great amount of filtered projection and image 
 data .
the cpu data access mechanism with multi - level caches is 
 illustrated in figure 7 .
specifically , a cache can be 
 used to reduce the average time to access data in the main memory 
 ( ram ) .
the cache is a smaller , faster memory chip which stores 
 copies of data from the most frequently used main memory 
 locations .
as long as a majority of memory accesses are to the 
 cached memory locations , the average latency of memory accesses 
 will be reduced to the cache latency , instead of the main memory 
 latency .
the l2 cache is faster than ram and about 
 1    2 mb .
the slowest ram is 1  4 gb . when the processor needs to read from or write to the main memory , 
 it first checks if the data is in the cache .
if it is in the 
 cache , we say that a cache hit has occurred ; otherwise a cache 
 miss is counted . in the case of a cache hit ,
however , in the case of a 
 cache miss , it takes much longer time to access the data .
due to 
 the limited cache capacity , one way to execute the code 
 efficiently is to increase the hit rate by optimizing the data 
 structures .
usually , projection data are sequentially stored in the order of 
  , while reconstructed images are stored rowwisely .
thus , 
 for implementation of the right - angle and symmetry , the access to 
 8 pairs of projection and image data will very likely result in 
 cache misses due to the address gaps , as shown in 
 figure 8 .
such misses within the inner loop will 
 cause a significant latency . to address this problem , in our 
 optimized data structures all the data are arranged into blocks 
 indexed to reflect symmetric relationships .
therefore , the cache 
 miss rate can be greatly reduced in the inner loop .
the simd technique enables the data - level parallelism like in a 
 vector processor , as shown in figure 9 . with an simd 
 processor ,
one instruction can process a block of data at a time 
 instead of just one datum , which is much more efficient than the 
 conventional single instruction single - data ( sisd ) technique . 

small - scale ( 64 or 128 bits ) simd operations are now popular 
 supported by general pc cpus , such as those from intel and amd 
 [ 12 , 13 ] .
we use the intel sse ( streaming simd extensions ) 
 instruction set to implement the simd technique in our 
 backprojection process . within the inner loop ,
we backproject 8 
 projection data onto 8 pixels , according to the same 
 instructions such as interpolation , weighting , and accumulation . 

therefore , we have a perfect situation to employ the simd 
 technique . as the sse only supports simultaneous processing of 4 
 floating data at a time ( 128-bits register ) , 8 data are processed 
 in two groups . in recent years
their efforts have now 
 shifted from improving the clock speed to increasing the number of 
 cores within a processor .
however , a processor with more 
 than one core can not achieve a better performance unless parallel 
 computation schemes are applied .
therefore , to take advantage of 
 multi - core processors , multi - threaded programing must be done . from the flowchart of our algorithm ,
thus , the backprojection can be 
 implemented in parallel by assigning different loop ranges to 
 various cores of the processor .
after all the threads are 
 finished , the final result can be assembled from the results of 
 each thread . in our implementation
( figure 4(b ) and 
 algorithm 2 ) , we divide the loop of 
 y , instead of the loop of . usually , the 
 number of cores on a pc is 2 , 4 or 8 , it is not common for 
 n/4 to be an integer , but it 
 is always the case for 
 ny to be 256 , 512 , 
 or 1024 .
our parallel implementation on a multi - core pc is more 
 efficient than that on a pc cluster in terms of time required for 
 data exchange between threads . in our case , the data exchange is 
 via on board ram bus , while the pc cluster 's data exchange 
 via local network is significantly slower .
the intel c++ compiler creates applications that can run at the 
 fastest speeds on the intel processors .
it can take the full 
 advantage of the intel processors when compiling codes and 
 generating object files .
this provides an integrated development 
 environment . in our implementation , we use it to optimize the code 
 for pentium d and xeon processors .
to test the gain of our scheme , we ran our accelerated code on 
 pentium d and xeon pcs .
, different sizes of projection 
 datasets and reconstructed images were tested to evaluate the 
 efficiency under various conditions . here to test effeteness of each technique , the reconstruction 
 times and speedups are tested by applying them gradually .
the 
 reconstruction experiments are done based on our hp6200 
 workstation and reconstructing a 512  512 image from a 
 projection dataset ( 1160  672 ) .
the reconstruction results 
 are shown by applying techniques step by step 
 ( table 2 ) .
the overall speedup and individual speedups 
 for each technique are also calculated to show the efficiency of 
 them .
the speedup results for different projection datasets and image 
 matrix sizes are shown in tables 3 , 4 , and 
 figure 10 .
significant speedups were achieved using our scheme . 
 in the case of 1160 views and 512  512 image , the 
 reconstruction time was decreased from 52 seconds to 1.35 seconds . 
 for a one - core computer ,
the speedup was more than 20 times . for a 
 two - core computer ,
the images reconstructed using our scheme and the 
 conventional method are shown in figure 11 .
all the 
 images were reconstructed using 32-bit floating - point data and 
 were displayed in the same window [ 0.97 , 1.05 ] . the images 
 reconstructed using our accelerated and conventional schemes are 
 essentially the same .
as the ct reconstruction algorithm is highly parallelizable , the 
 speedup can be improved with more cores almost linearly . for 
 example , with two quad - core processors , the speedup that could be 
 achieved is more than 100 .
as compared to other acceleration 
 techniques , such as those based on specialized hardware and 
 graphics cards , our general purpose pc - based scheme is much 
 cheaper without compromising image quality .
for example , a 
 general purpose hp or dell workstation with a top - line two 
 quad - core processor and 8 gb ram is less than $ 7000 .
all 
 calculations are based on 32-bit floating point data , providing 
 sufficient accuracy for medical imaging applications . in terms of the absolute reconstruction time for a 512  512 
 image from 1160 projection views
if the latest multi - core processor is 
 used , the total time can be easily decreased by several folds .
as 
 the computers we have still use the previous generation processor , 
 the potential improvement is at least 5 times if we are equipped 
 with the latest quad - core processors , that is , the 
 reconstruction time may be reduced to 0.3 second .
hence , it is 
 quite promising for real - time ct applications , such as project on 
 bolus - chasing ct angiography . in conclusion ,
our acceleration scheme has integrated several 
 techniques including utilization of geometric symmetry , 
 optimization of data structures , single - instruction multiple - data 
 ( simd ) processing , multi - threaded computation , and an intel c++ 
 complier . as a result
, it has speeded up the reconstruction 
 process by 40 times , as compared to the conventional 
 implementation on a general purpose pc with 2 cores .
further work 
 is in progress to improve our results using the latest pcs and 
 extend our scheme for cone - beam reconstruction .
 expensive computational cost is a severe limitation in ct reconstruction for clinical applications that need real - time feedback . 
 a primary example is bolus - chasing computed tomography ( ct ) angiography ( bca ) that we have been developing for the past several years . to accelerate the reconstruction process using the filtered backprojection ( fbp ) method , specialized hardware or graphics cards can be used . however , specialized hardware is expensive and not flexible . 
 the graphics processing unit ( gpu ) in a current graphic card can only reconstruct images in a reduced precision and is not easy to program . in this paper , 
 an acceleration scheme is proposed based on a multi - core pc . in the proposed scheme , 
 several techniques are integrated , including utilization of geometric symmetry , optimization of data structures , single - instruction multiple - data ( simd ) processing , multithreaded computation , and an intel c++ compilier . 
 our scheme maintains the original precision and involves no data exchange between the gpu and cpu . 
 the merits of our scheme are demonstrated in numerical experiments against the traditional implementation . 
 our scheme achieves a speedup of about 40 , which can be further improved by several folds using the latest quad - core   processors . 