image quality greatly depends on how accurately physical phenomena are modelled in the system , which is represented as the system response matrix ( srm ) [ 1 , 2 ] .
the srm models the probability of detection , of an annihilation produced in voxel b , in a detector element or crystal pair d , pd , b .
the size of the srm is imposed by the number of voxels n that comprises the field of view ( fov ) and by the number of detector elements of the scanner , j. the size of the srm , n  j , is significant due to the large number of crystals used by high resolution scanners and the fine granularity of the fov .
fortunately , most elements pd , b in an srm are typically zero , which means that the number of nonzero srm elements to store is significantly lower than n  j. the srm can be directly obtained by taking measurements where a point source is placed in different positions over the fov , interpolating intermediate positions when needed .
this produces a highly accurate srm at the cost of long acquisition times [ 4 , 5 ] . to greatly reduce acquisition time
however , rather than directly measuring the system matrix such methods are used only to estimate the shift - variant point spread function .
the speed of these methods outperforms others , at the expense of limited precision [ 7 , 8 ] .
these methods produce a more accurate srm than analytical methods , but do not include the physical phenomena that are only found by measuring the srm . regarding the time necessary to calculate the srm , a simulation is longer than analytical computations , but less time consuming than that necessary to obtain measurements experimentally .
the high granularity of current pet scanners , used to obtain the best possible spatial resolution , makes the calculation of the srm elements a cumbersome process .
the use of cylindrical symmetries , taking advantage of the polygonal architecture of pet scanners , is a common approach to reduce the number of srm elements required [ 915 ] . in the specific case of a monte carlo - based srm ,
although the most common model used to represent the image space is the cubic voxel , there exists a number of alternative basis functions .
the use of polar voxels provides a convenient model to exploit scanner symmetries , which significantly reduces the simulation time and space necessary to store the srm only modelling a small portion of the fov .
alternatively , spherical basis functions ( blobs ) have been shown to represent a more suitable basis to model the continuous radiotracer distribution . improved
noise performance has been demonstrated , compared to cubic and polar voxels [ 1820 ] .
moreover , better spatial resolution can be obtained with spherical basis functions , compared to postreconstruction filtered cubic and polar voxels .
therefore , the number of basis functions intersected by any given line of response ( lor ) is higher using blobs compared to voxels .
this characteristic produces that the srm has a significantly large number of nonzero srm elements , resulting in long cpu - based reconstruction times .
the forward and backward projections in a maximum - likelihood- ( ml- ) based reconstruction algorithm perform independent operations on each lor and voxel , respectively , which makes these operations highly parallelizable . graphical processing units ( gpus ) technology have been successfully employed to enhance speed of image reconstruction in both precomputed [ 2123 ] and on - the - fly [ 24 , 25 ] srms , which makes it a potential candidate for blob - based reconstruction as well .
an alternative approach to parallelize the reconstruction process is to use cpu - multicore architectures . in this work
we have implemented an optimized cpu - multicore version of the reconstruction algorithm using the message - passing - interface ( mpi ) libraries ( http://www.mcs.anl.gov/research/projects/mpi/index.htm ) .
this paper explores the suitability of these two approaches for the special case of ml - expectation - maximization ( em ) pet reconstruction using a blob - based monte carlo precomputed srm .
the use of a large precomputed srm in gpu technology represents one of the major challenges of this work .
in this study the small animal scanner madpet - ii , shown in figure 1 , is used as a model for all the simulations carried out .
madpet - ii has a radial diameter of 71  mm and an axial fov of 18.1  mm .
the ring contains 18 modules , where each module has two layers of 8  4 lyso crystals with individually read out electronics based on avalanche photodiodes .
the size of the crystals in the front layer is 2  2  6  mm and that of the rear layer is 2  2  8  mm .
the dual layer provides information of the depth of interaction , which mitigates parallax error . the graphics card used in image reconstruction was the nvidia tesla c2070 , based on a fermi architecture .
the gpu contains 448 cores ( thread processors or shaders ) running at 1.15  ghz ( 575  mhz core ) , 6  gb gddr5 on - board global memory , 48  kb of shared memory per block , and a bandwidth of 144  gb / s .
threads are internally organized as a grid of thread blocks , where threads in the same block communicate through shared memory .
the maximum number of blocks is 65  535 and the maximum number of threads per block is 1 024 .
the card has a peak performance of 515.2  gflops / s for double precision operations .
the algorithm was implemented for both the gpu and an eight - core desktop pc ( intel(r ) core(tm ) i7 cpu 950 @ 3.07  ghz ) with 12  gb of ram .
an optimized algorithm allowed the use of either a single core or multiple cores and distributed memory .
the total number of crystals of madpet - ii is 1  152 ; for image reconstruction using cubic voxels an image space of 140  140  40 ( voxels of size 0.5  0.5  0.5  mm ) was used , resulting in j  n = 10 elements .
however , the srm is highly sparse , significantly reducing the number of elements necessary for storage .
a popular method to reduce the number of srm elements to calculate and to reduce the srm storage size is to exploit the scanner symmetries .
the typical polygonal architecture in pet scanners allows the use of polar symmetries . for a monte carlo - based srm
this approach allows not only the reduction of the file size of the srm , but also the simulation time .
the number of symmetries in a pet system is intrinsically linked to its geometry , specifically to the number of detector modules employed over the circumference .
the high resolution small animal scanner madpet - ii contains 18 block detectors . using 18 rotational symmetries and 2 reflection symmetries transaxially and 2 reflection symmetries axially , up to 72 symmetries can be exploited .  
therefore , a factor of 72 reduction in the simulation time and in the srm file size is achieved .
more details about this implementation of symmetries in madpet - ii can be found in .
the srm was calculated using gate   in the local grid facility which had 108 nodes ( two quad core xeon e5420 @ 2.50 ghz machines ) for sequential processing .
the total number of simulated events in the srm was 3.7  10 , which represents 4.3  10 events per voxel .
the simulation of the 3.7  10 events was split in 200 parallel jobs , where each parallel job took approximately 12 hours .
the srm was simulated using a back - to - back gamma source , thus half life , noncollinearity and positron range were ignored .
singles list mode data were processed after simulation to select coincidences using a coincidence timing window of 20  ns .
an estimation of the continuous radiotracer distribution f(r ) is represented in image space as a linear combination of image coefficients and basis functions , expressed as 
 ( 1)f(r)f^(r)=i=1nci(rri ) , 
 where r represents the space coordinate , f^(r ) the estimated distribution , ci the image coefficients , (r ) the basis functions , n the total number of voxels , and ri the placement grid .
spherical basis functions have compact support , that is , the function is zero beyond a given value ( blob radius ) , but have a smoother behaviour than the traditional cubic voxels .
the spherical basis functions used in this work are based on the kaisser - bessel function , described as 
 ( 2)(r)m , a,=1im()[1(ra)2]mim(1(ra)2 ) , 
 where m is the order of the bessel function i , a is the blob radius , and  is the taper parameter .
a thorough investigation of these parameters is found in   in which the optimal values are obtained .
the study , performed in frequency space , determines that the optimal value of m is 2 ( the first derivative is continuous at the blob boundary ) , the optimal radius a is 1.994 (  being the distance between elements of the underlying grid ) , and the optimal  is 10.4 . 
 figure 2 shows the wedge - like source simulated to calculate the srm elements using symmetries , with an underlying polar grid and blobs placed over the grid using a body - centred strategy with the polar voxels used as reference . using spherically symmetric basis functions and exploiting spherical symmetries ,
the final srm was 5.3  gb in size and contained 7.04  10 nonzero elements ( out of the 5.79  10 : 8740 blobs   662  976 lors ) .
if eight cubic symmetries were used ( 2 rotational + 2 reflectional transaxially + 2 reflectional axially ) , instead of spherical symmetries , the resulting srm file size would be 47.2  gb with the same statistical quality as that used .
all phantoms were simulated using a back - to - back gamma source in agreement with the simulation of the srm .
simulated data is stored in lor - histograms where each histogram bin corresponds to an lor , and no preprocessing was applied to the data [ 9 , 29 ] .
a homogeneously filled phantom with a hot and a cold rod inserts has been simulated to study the image quality .
the cylinder is 20  mm long and 30  mm radius , while the rod inserts are 20  mm long and 10  mm radius .
the simulated activity of the phantom was 32  mbq ( 0.86  mci ) and 74  mbq ( 2  mci ) in the hot and warm regions , respectively ( 3  :  1 ratio ) .
the experiment time simulated was 155 seconds , producing a total of 4.5 10 coincidences . to study the spatial resolution , an ellipsoidal phantom with six hot point sources and six cold spheres ( 1.5  mm  ) , placed radially along the ellipsoid and separated by 5  mm on a warm background , was simulated .
the ellipsoid was 35  mm by 20  mm transaxially and 2  mm long .
the phantom was placed 12.5  mm off - centre covering more than one half of the fov .
the activity simulated in the point sources was 3.7  mbq ( 0.1  mci ) and 207.2  mbq ( 5.6  mci ) in the background ( 200  :  1 ratio ) , producing a total of 9.6 10 coincidences .
the background was used in order to mitigate the resolution enhancement caused by the nonnegativity constraint of ml - em [ 5 , 30 ] .
the digital mouse phantom moby   was simulated with a total activity of 0.2  mci .
five bed positions were necessary to acquire the whole mouse with six overlapping slices ( 3  mm ) between bed positions , with a time scan of four minutes per bed position .
the activity simulated in each organ is the default relative activity set by the moby phantom files .
it was not a purpose of this study to simulate realistic activity concentrations in each organ .
the noise performance obtained with each of the implementations explored was measured in the image quality phantom using three figures of merit , the coefficient of variation ( cv ) , the contrast to noise ratio ( cnr ) , and the correlation coefficient ( cc ) .
the spatial resolution was also studied by measuring the full width at half maximum ( fwhm ) across the point sources of the ellipsoidal phantom described in section 2.5 .
these figures of merit were measured to compare the different implementations presented here , and not to assess the performance of spherical basis functions .
the cv and cnr were measured over regions of interest ( roi ) of size 8  9  0.5  mm , placed far enough from the boundaries so that there were no edge effects .
the cv is commonly used as a normalized measurement of noise in a given roi and is described as 
 ( 3)cv=ss , 
 where s is the mean value and s is the standard deviation measured in the roi .
the cnr is a measure of noise performance between two rois   given by 
 ( 4)cnr=sb(s2+b2)/2 , 
 where b is the mean value and b is the standard deviation measured in the background .
the cc between two images is a statistical similarity measure defined as 
 ( 5)ccab=i=1n(f^aifa)(f^bifb)(i=1n(f^aifa)2)(i=1n(f^bifb)2 ) , 
 where f^ai is the intensity value at voxel i in image a , fa is the mean value of image a , and similarly for image b. a cc value of 1 represents two perfectly correlated images , while a value of 1 represents two completely uncorrelated images .
the spatial resolution was measured as the fwhm taken from a profile drawn across the hot point sources embedded in the ellipsoidal phantom .
one of the main challenges addressed in this work was the difficulty to cope with the amount of information needed by a single thread to perform both forward and backward projection , given that the srm was too large to be stored in local memory , registers or shared memory . while a memory access to a register or shared memory takes only one clock cycle , access to global memory takes 400600 clock cycles and access to texture memory takes 1100 clock cycles ( depending on cache locality ) .
the floating point values used in this work were distributed as follows : given the significant size of the srm ( 5.3  gb ) this was stored in global memory , while the arrays stored in texture memory were the ratio between measurements and projected image estimate ( 2.52  mb ) , a look - up table used to unfold the symmetries ( 2.40  mb ) and the image estimate used in the forward projection ( 2.40  mb ) .
the number of blocks in the gpu grid was optimized to achieve the shortest reconstruction time per iteration , being 32 blocks obtained empirically .
therefore , memory access to the srm elements in the forward projection operation was consecutive , as opposed to the backward projection , where access to the srm elements was highly irregular .
the srm organization implies that the forward projection was a gathering operation , while the backward projection was a scattering operation , which is slower . during the forward projection ,
each thread reads from a number of blob coordinates and calculates the projection of one lor , hence writing in different memory positions .
however , during the backward projection , each thread reads from a number of measurements , back - projects the current estimate to image space , and writes in the corresponding blob coordinate .
given that threads are organized by lors , different threads can write in the same memory position .
this represents what is known as a race condition . to avoid such problem three different approaches
an srm ordered by consecutive lor indices ( detection ) was used in the forward projection , while an srm ordered by consecutive blob indices ( image ) was used in the backward projection .
this second srm represents the transpose of the initial srm although stored in a different file due to the use of sparse format . using this approach
the advantages are that the access to srm elements was in consecutive memory positions for both operations and that there was no overwrite hazards .
the disadvantage of this implementation is the time required for each iteration to load and unload each srm file in the gpu . to refer to this implementation
an srm ordered by consecutive lor indices was used in both , the forward and backward projections .
this approach represented an overwrite hazard during the backward projection when lors intersecting the same voxels were processed by threads in parallel .
this situation is likely to happen in those areas where most lors intersect , that is , the centre of the fov .
atomic writing prevents this situation from occurring , at the cost of longer computation time . to refer to this implementation
similarly to the implementation above , an srm ordered by consecutive lor indices was used in both , the forward and backward projections . to mitigate the speed problem caused by atomic operations ,
a combination of two strategies was followed : 
 ( a)the lors were sent to the threads using a random ordering , hence not following any spatial correlation . by introducing spatial randomness in the execution of lors in parallel ,
however , writing in the same memory positions still happened in a region at the centre of the fov;(b)atomic writing was exclusively used for those voxels located in the two central slices . otherwise a nonatomic operation was performed . 
 
the lors were sent to the threads using a random ordering , hence not following any spatial correlation . by introducing spatial randomness in the execution of lors in parallel ,
writing in the same memory positions still happened in a region at the centre of the fov ; atomic writing was exclusively used for those voxels located in the two central slices . otherwise a nonatomic operation was performed . 
 to refer to this implementation the term srmrand - at was used .
this approach used distributed memory , hence a portion of the srm is sent to each core to perform the forward and backward projections , respectively .
the srm ordered by lor indices was used in the forward projection , while a blobs - ordered srm is used for the backward projection , similar to the gpu implementation using srmreload .
subsequently , the srm was loaded and unloaded each time a forward and a backward operation was performed , adding a computational overhead .
the reconstruction time obtained in this work using gpu technology and spherical basis functions is comparable to the reconstruction time obtained using polar or cubic basis functions on a cpu .
however , the time performance is not as high as that published in other works [ 22 , 25 ] due to the nature of this approach , that is , the use of a large monte carlo precomputed srm .
this implies that the number of global memory accesses by each thread in a forward / backward projection corresponds to the number of nonzero elements of each lor / blob multiplied by the number of symmetries .
table 1 shows a comparison between the time performance measured using a single intel(r ) core(tm ) i7 cpu 950 @ 3.07  ghz , that is , nonparallelized ( np ) , with the time performance measured using the cpu - multicore implementation for 1 ( mpi-1 ) , 2 ( mpi-2 ) , 4 ( mpi-4 ) , and 8 ( mpi-8 ) cores , and finally the time performance measured with the gpu , using the different implementations described above ( srm ignoring the overwrite hazard , srmreload , srmatomic , and srmrand - at ) .
the results shown in table 1 were obtained using the image quality phantom described in section 2.5 . nevertheless , small variations were observed between different phantoms .
the cpu - single - core implementation has been taken as reference to calculate the improvement factors of the parallelized implementations .
table 1 shows that the gpu - based implementation where the overwrite hazard is ignored is the fastest implementation because atomic writing is not used , at the cost of producing unacceptable artefacts in the final reconstructed images ( see figure 3(c ) ) . using srmreload , memory access to srm elements
is performed consecutively , both in the forward and backward projection . in every other implementation ,
consecutive memory access is performed only in the forward projection but not in the backward projection , hence decreasing the time performance . using the implementation with srmreload the forward projection takes 97  s and 58  s to load the srm , while the backward projection requires only 87  s and 57  s for loading .
this represents an iteration time of 184  s if we consider only processing time . however , an extra 115  s is taken to load the srm in each operation .
however , by using the atomic operation only for those critical voxels where the probability of over - writing is high , the backward projection time of implementation using srmrand - at is reduced to 132  s , close to the backward projection time measured in the implementation where the overwrite hazard is ignored ( 121  s ) . moreover
, the artefacts obtained in the reconstructed images when the overwrite hazard is ignored are removed using srmrand - at ( figure 3(d ) ) . when enhancing speed of an image reconstruction algorithm , it is of critical importance to produce the same image for each implementation . to demonstrate that the implementations detailed in this work do not have an impact on image quality , the phantoms described in section 2.5 have been reconstructed using the cpu - single - core implementation , the cpu - multicore implementation using eight cores , and the gpu - based implementation using srmrand - at listed in table 1 .
the impact on the noise performance has been assessed using the image quality phantom described in section 2.5 , which has been reconstructed after 300 iterations using four of the implementations studied ( figure 3 ) : the cpu - single - core , the cpu - multicore , and two of the gpu - based implementations , without atomic operation and with atomic operation used only in the two central slices ( srmrand - at ) .
special mention is required for the gpu - based implementation without atomic writing ( figure 3(c ) ) where significant artefacts are observed , mainly in the centre of the fov .
as explained above , these artefacts are due to parallel threads overwriting in the same memory positions during the backward projection , due to the high overlapping between lors in the centre of the fov .
the images obtained with srmreload and srmatomic ( not shown in this work ) are practically identical to the one obtained with srmrand - at .
a profile across the four reconstructed phantoms is shown in figure 3(e ) , demonstrating great resemblance between the images obtained with the cpu - single - core , the cpu - multicore , and the gpu - based implementation with atomic writing , while the artefact observed in figure 3(c ) is clearly visible in figure 3(e ) in the black profile .
for quantitative assessment , the cv ( figure 4 ) was measured in the hot , warm , and cold rois for the three different implementations , every 10 iterations for 300 iterations .
the cv at iteration 300 is 0.15 and 0.28 for the hot and warm rois , respectively , for all three implementations , while the cv in the cold roi for the cpu - based implementations is 0.72 and for the gpu - based implementation 0.73 . in all cases
the cv follows an increasing trend due to the known noise increase as more ml - em iterations are calculated .
differences between the cpu - based implementations ( single - core and multicore ) are below 0.08% ( figure 4(b ) ) , while higher differences were observed between the cpu - single - core and gpu - based implementations , where differences between 0.37% ( hot roi ) and 0.81% ( cold roi ) were measured .
these differences are due to the different floating point precisions available in the cpu and the gpu .
while the precision had little effect on individual srm elements , the cumulative effect produced differences in the reconstructed images .
however , these differences are expected to be dominated by statistical errors in the data . 

figure 5 shows the evolution of the cnr measured between the hot and warm rois for 300 iterations at every 10 iteration .
the differences between the cpu - based and the gpu - based implementations are shown in figure 5(b ) . similarly to the cv study , differences between the cpu - based implementations ( single - core and multicore ) are below 0.02% , while small differences are observed between the cpu - single - core and the gpu implementations .
finally , the cc ( figure 6 ) was measured between the entire volumes of the resulting reconstructed image quality phantoms obtained with the cpu - based implementation , the mpi - based implementation , and the gpu - based implementation .
similarly to the cv and the cnr studies , the cc measured between the cpu - based and the mpi - based implementations shows perfect correlation , while the comparison between the gpu - based implementation with the cpu - based and the mpi - based implementations show high correlation initially , but the trend is to slightly decrease at later iterations .
however , the cc at 300 iterations is over 98.5% so that the difference can be considered negligible . to confirm that the trend was not exacerbated as more iterations are computed ,
the cc between the mpi - based implementation and the gpu - based implementation was computed over 600 iterations , showing a cc of 94% at 600 iterations and a slight trend correction . 

figure 7 shows the spatial resolution phantom reconstructed after 300 iterations using the cpu - single - core , the cpu - multicore , and the gpu - based implementations ( using srmrand - at ) .
the profile drawn across the point sources in the three reconstructed phantoms , shown in figure 8 , confirms the conclusions made in the noise study .
the cpu - single - core and the cpu - multicore implementations produce identical results while the images reconstructed using the gpu implementation are slightly different .
this is further confirmed by table 2 , where the fwhm ( mm ) measured from each point source and each profile is shown ( ps 1 corresponds to the point source located closer to the centre of the fov while ps 6 corresponds to the point source located closer to the edge of the fov ) .
the spatial resolution decreases as the point source is located farther from the centre of the fov due to the parallax effect .
similarly to the noise assessment study , the fwhm was measured over 300 iterations at steps of 10 iterations ( figure 9(a ) ) , showing small differences measured between the cpu - based and the gpu - based implementations below 0.4% at 300 iterations ( figure 9(b ) ) . for qualitative assessment , the moby
phantom   has been reconstructed using blobs in the gpu and using the cpu with the mpi libraries . for comparison purposes ,
however , it is important to highlight that the focus of this work is not to compare these two basis functions .
spherical basis functions provide better noise performance than cubic voxels , so in order to perform a fair comparison , the image reconstructed using cubic voxels was filtered to match the noise performance of that achieved using spherical basis functions , which produces a visible image detail degradation . 

figure 10(b ) shows the reconstructed phantom obtained using cubic voxels with a postreconstruction gaussian filter of  = 0.5  mm for comparison purposes .
a  = 0.5  mm for the gaussian filter has been applied to match the noise levels obtained with cubic voxels and blobs .  
figure 10(c ) shows the phantom reconstructed with spherically symmetric basis functions using the gpu implementation and figure 10(d ) shows the phantom reconstructed with spherically symmetric basis functions using the cpu - multicore implementation .
300 iterations were used to reconstruct each bed position in the three reconstructed moby phantoms presented here .
it can be noticed that the thyroids and brain are more visible , and boundaries better delineated using spherical basis functions , compared to filtered cubic voxels ( figure 10(b ) ) , as shown in figures 10(c ) and 10(d ) .
the reconstruction time necessary using cubic voxels in a single core in the cpu ( figure 10(b ) ) was 110 hours , the cpu - multicore implementation using spherical basis functions ( figure 10(d ) ) was 220 hours , and the gpu - based implementation using spherical basis functions ( figure 10(c ) ) was 88 hours .
if the same phantom was reconstructed using spherical basis functions on a single core , the necessary reconstruction time would be 380 hours . from figure 10
it can be seen that the combination of blob - based reconstruction and a precomputed monte carlo srm using gpu technology is a feasible alternative , not only for simple phantom geometries as those shown in section 4.2 , but also for multistage and highly realistic phantoms as the moby phantom .
increasing granularity in pet scanners provides improved spatial resolution while increasing the number of detector elements at the same granularity improves sensitivity .
however , increased resolution means that the calculation of the srm is an extremely cumbersome task , particularly for simulation - based srm calculation . nevertheless , monte carlo - based system matrices for iterative statistical image reconstruction applied to emission tomography are growing in popularity due to their image quality advantages .
the extended availability of affordable computing power means that significant efforts are being put into sophisticated improvements of the system response model .
overlapping spherically symmetrical basis functions have clear advantages over nonoverlapping ( cubic or polar voxels ) even at the cost of a significantly high number of nonzero elements in the srm , resulting in large srm file sizes and long reconstruction times . these problems can be partly overcome by exploiting cylindrical symmetries to reduce the simulation time , the number of nonzero srm elements , and hence the file size necessary to store the srm .
the combination of spherically symmetric basis functions and cylindrical symmetries makes this approach feasible for use in a clinical or preclinical application . however , reconstruction time is then the main concern due to the still large number of nonzero srm elements required to process the forward and backward projection .
ordered - subsets- ( os- ) em represents a common approach to speed up the reconstruction process .
while os - em can be implemented in gpu technology , it requires a device - dependent level of complexity , and its inclusion may reduce the generality of the study presented here .
moreover , subset choice interacts with both speed and image quality while hardware - based solutions decouple this relationship .
while it is expected that the combination of os - em and gpu technology can effectively further reduce reconstruction times , this is beyond the scope of this investigation .
this work presents a srm - generic hardware implementation that achieved reconstructed images 4.3 times faster using gpu technology compared to an optimized cpu - single - core implementation and 2.5 times faster than a cpu - multicore ( 8) implementation .
a cpu - multicore implementation decreased the reconstruction time by only 1.7 times compared to the single - core implementation .
differences in image performance were assessed from an image quality and a spatial resolution perspective .
negligible differences were demonstrated between cpu - based single - core and multicore implementations , and small differences between cpu - single - core and gpu - based implementations were observed .
differences below 1% for the cv , below 2.5% for the cnr , and below 0.4% for the spatial resolution were measured between the cpu - single - core implementation and the gpu - based implementation .
 statistical iterative methods are a widely used method of image reconstruction in emission tomography . 
 traditionally , the image space is modelled as a combination of cubic voxels as a matter of simplicity . after reconstruction 
 , images are routinely filtered to reduce statistical noise at the cost of spatial resolution degradation . 
 an alternative to produce lower noise during reconstruction is to model the image space with spherical basis functions . 
 these basis functions overlap in space producing a significantly large number of non - zero elements in the system response matrix ( srm ) to store , which additionally leads to long reconstruction times . 
 these two problems are partly overcome by exploiting spherical symmetries , although computation time is still slower compared to non - overlapping basis functions . in this work , 
 we have implemented the reconstruction algorithm using graphical processing unit ( gpu ) technology for speed and a precomputed monte - carlo - calculated srm for accuracy . 
 the reconstruction time achieved using spherical basis functions on a gpu was 4.3 times faster than the central processing unit ( cpu ) and 2.5 times faster than a cpu - multi - core parallel implementation using eight cores . 
 overwriting hazards are minimized by combining a random line of response ordering and constrained atomic writing . 
 small differences in image quality were observed between implementations . 
1. Introduction
2. Materials and Methods
3. Hardware Implementations
4. Results and Discussion
5. Conclusions