image denoising is commonly used in medical imaging in order to help medical doctors to see abnormalities in the images .
image denoising was first applied to 2d images [ 13 ] and then extended to 3d data [ 46 ] , 3d data can either be collected as several 2d images over time or as one 3d volume .
a number of medical imaging modalities ( e.g. , computed tomography ( ct ) , ultrasound ( us ) and magnetic resonance imaging ( mri ) ) now provide the possibility to collect 4d data , that is , time - resolved volume data .
this makes it possible to , for example , examine what parts of the brain that are active during a certain task ( functional magnetic resonance imaging ( fmri ) ) . while 4d ct data makes it possible to see the heart beat in 3d
, the drawback is that a lower amount of x - ray exposure has to be used for 4d ct data collection , compared to 3d ct data collection , in order to not harm the patient . when the amount of exposure is decreased , the amount of noise in the data increases significantly .
three - dimensional image denoising has previously been applied to several time points independently , but there has not been much work done on true 4d image denoising where the algorithm considers several volumes at the same time ( and not a single volume at a time ) .
applied 4d anisotropic diffusion filtering to ultrasound volumes and jahanian et al .   applied 4d wavelet denoising to diffusion tensor mri data . for ct data , it can be extra beneficial to use the time dimension in the denoising , as some of the reconstruction artefacts vary with time .
it is thereby possible to remove these artefacts by taking full advantage of the 4d data .
while true 4d image denoising is very powerful , the drawback is that the processing time increases exponentially with respect to dimensionality .
the rapid development of graphics processing units ( gpus ) has resulted in that many algorithms in the medical imaging domain have been implemented on the gpu , in order to save time and to be able to apply more advanced analysis . to give an example of the rapid gpu development ,
. some examples of fields in medical imaging that have taken advantage of the computational power of the gpu are image registration [ 913 ] , image segmentation [ 1416 ] and fmri analysis [ 1720 ] . in the area of image denoising ,
howison   made a comparison between different gpu implementations of anisotropic diffusion and bilateral filtering for 3d data .
su and xu   in 2010 proposed how to accelerate wavelet - based image denoising by using the gpu .
describe gpu - based image manipulation and enhancement techniques for dynamic volumetric medical image visualization , but enhancement in this case refers to enhancement of the visualization , and not of the 4d data .
recently , the gpu has been used for real - time image denoising . in 2007 , chen et al .   used bilateral filtering   on the gpu for real - time edge - aware image processing .
fontes et al .   in 2011 used the gpu for real - time denoising of ultrasound data and goossens et al .   in 2010 managed to run the commonly used nonlocal means algorithm   in real time .
to our knowledge , there has not been any work done about true 4d image denoising on the gpu . in this work we therefore present a novel algorithm , based on local adaptive filtering , for 4d denoising and
describe how to implement it on the gpu , in order to decrease the processing time and thereby significantly increase the clinical value .
   in this section , the algorithm that is used for true 4d image denoising will be described . to show how a higher number of dimensions , the power of dimensionality , can improve the denoising result ,
the size of the 4d data is 127  127  9  9 , but there is no signal variation in the last two dimensions .
the data contains a large step , a thin line , and a shading from the top left corner to the bottom right corner .
the denoising was done on one 127  127 image , for the 3d case , the denoising was done on one 127  127  9 volume and for the 4d case all the data was used .
a single anisotropic lowpass filter was used for the denoising , and the filter had the same dimensionality as the data and was oriented along the structures .
the original test data , the test data with noise and the denoising results are given in figure 1 .
it was introduced for 2d by knutsson et al . in 1983   and then extended to 3d in 1992 . in this work ,
the main idea is to first estimate the local structure tensor   ( by using a first set of filters ) in each neighbourhood of the data and then let the tensor control the reconstruction filters ( a second set of filters ) .
the term reconstruction should in this paper not be confused with the reconstruction of the ct data .
the local structure tensor t is in 4d a 4  4 symmetric matrix in each time voxel , 
 ( 1)t=(t1t2t3t4t2t5t6t7t3t6t8t9t4t7t9t10)=(xxxyxzxtxyyyyzytxzyzzzztxtytzttt ) , 
							 and contains information about the local structure in the data , that can be used to control the weights of the reconstruction filters .
the result of the adaptive filtering is that smoothing is done along structures ( such as lines and edges in 2d ) , but not perpendicular to them .
compared to more recently developed methods for image denoising ( e.g. , nonlocal means , anisotropic diffusion   and bilateral filtering ) , adaptive filtering is in our case used for 4d image denoising for three reasons .
nonlocal means can give very good results , but the algorithm can be extremely time consuming ( even if gpus are used ) .
bilateral filtering does not only require a multiplication for each filter coefficient and each data value , but also an evaluation of the intensity range function ( e.g. ,   an exponential ) which is much more expensive to perform than a multiplication .
second , the tuning of the parameters is for our denoising algorithm rather easy to understand and to explore .
when a first denoising result has been obtained , it is often obvious how to change the parameters to improve the result .
third , the adaptive filtering approach has been proven to be very robust ( it is extremely seldom that a strange result is obtained ) .
adaptive filtering has been used for 2d image denoising in commercial clinical software for over 20 years and a recent 3d study   proves its potential , robustness , and clinical acceptance .
the nonlocal means algorithm only works if the data contains several neighbourhoods with similar properties .
the local structure tensor can , for example , be estimated by using quadrature filters [ 5 , 30 ] .
quadrature filters q are zero in one half of the frequency domain ( defined by the direction of the filter ) and can be expressed as two polar separable functions , one radial function r and one directional function d , 
 ( 2)q(u)=r(||u||)d(u ) , 
							 where u is the frequency variable .
the radial function is a lognormal function 
 ( 3)r(||u||)=exp  ( c  ln  2(||u||u0 ) ) ,  c=4b2ln  ( 2 ) , 
							 where u0 is the centre frequency of the filter and b is the bandwidth ( in octaves ) .
the directional function depends on the angle  between the filter direction vector n^ and the normalized frequency coordinate vector u as cos  (  ) , 
 ( 4)d(u)={(utn^)2,utn^>0,0,otherwise . 
							
quadrature filters are cartesian nonseparable and complex valued in the spatial domain , the real part is even and in 2d acts as a line detector , while the imaginary part is odd and in 2d acts as an edge detector . in 3d , the even and odd filters correspond to a plane detector and a 3d edge detector . in 4d ,
the complex - valued filter response q is an estimate of a bandpass filtered version of the analytical signal with magnitude a and phase  , 
 ( 5)q = a(cos  ( )+isin  ( ))=aei. 
							 the tensor is calculated by multiplying the magnitude of the quadrature filter response qk with the outer product of the filter direction vector n^k and then summing the result over all filters k , 
 ( 6)t=k=1nf|qk|(c1n^kn^ktc2i ) , 
							 where c1 and c2 are scalar constants that depend on the dimensionality of the data [ 5 , 30 ] , nf is the number of quadrature filters and i is the identity matrix .
the resulting tensor is phase invariant , as the magnitude of the quadrature filter response is invariant to the type of local neighbourhood ( e.g. , in 2d bright lines , dark lines , dark to bright edges , etc . ) .
this is in contrast to when the local structure tensor is estimated by using gradient operators , such as sobel filters .
the number of filters that are required to estimate the tensor depends on the dimensionality of the data and is given by the number of independent components of the symmetric local structure tensor .
the required number of filters is thus 3 for 2d , 6 for 3d and 10 for 4d . the given tensor formula ,
it is possible to spread 6 filters evenly in 3d , but it is not possible to spread 10 filters evenly in 4d .
for this reason , 12 quadrature filters have to be used in 4d ( i.e. , a total of 24 filters in the spatial domain , 12 real valued and 12 complex valued ) . to apply 24 nonseparable filters to a 4d dataset
requires a huge number of multiplications . in this paper a new type of filters , monomial filters ,
monomial filters also have one radial function r and one directional function d. the directional part of the monomial filters are products of positive integer powers of the components of the frequency variable u. the monomial filter matrices of order one , f1 , and two , f2 , are in the frequency domain defined as 
 ( 7)f1,n = r(||u||)u^n ,  f2,mn = r(||u||)u^mu^n . 
							
the monomial filters are first described for 2d and then generalized to 4d . in 2d ,
the frequency variable is in this work defined as u = [ u  v ] .
the directional part of first - order monomial filters are x , y in the spatial domain and u , v in the frequency domain .
two - dimensional monomial filters of the first - order are given in figure 2 .
the directional part of second - order monomial filters are xx , xy , yy in the spatial domain and uu , uv , vv in the frequency domain .
the monomial filter response matrices q are either calculated by convolution in the spatial domain or by multiplication in the frequency domain . for a simple signal with phase  ( e.g. , s(x ) = acos  ( ux +  ) ) ;
the monomial filter response matrices of order one and two can be written as 
 ( 8)q1=iasin  ( )[uv]t , q2=acos  ( )(uuuvuvvv ) . 
									
the first - order products are odd functions and are thereby related to the odd sine function , the second order products are even functions and are thereby related to the even cosine function ( note the resemblance with quadrature filters that have one even real part and one odd imaginary part ) . by using the fact that u + v = 1 , the outer products of the filter response matrices give 
 ( 9)q1q1t = sin  2()|a|2(uuuvuvvv),q2q2t = cos  2()|a|2(uuuvuvvv ) . 
									 the local structure tensor t
, it is clear that the estimated tensor , as previously , is phase invariant as the square of one odd part and the square of one even part are combined . for information about how to calculate
a total of 14 nonseparable 4d monomial filters ( 4 odd of the first - order ( x , y , z , t ) and 10 even of the second - order ( xx , xy , xz , xt , yy , yz , yt , zz , zt , tt ) ) with a spatial support of 7  7  7  7 time voxels are applied to the ct volumes .
the filters have a lognormal radial function with centre frequency 3/5 and a bandwidth of 2.5 octaves .
the filter kernels were optimized with respect to ideal frequency response , spatial locality , and expected signal - to - noise ratio [ 5 , 33 ] . by using equation ( 10 ) for the 4d case , and replacing the frequency variables with the monomial filter responses , the 10 components of the structure tensor
are calculated according to 
 ( 11)t1=fr1fr1+fr5fr5+fr6fr6+fr7fr7  + fr8fr8,t2=fr1fr2+fr5fr6+fr6fr9+fr7fr10  + fr8fr11,t3=fr1fr3+fr5fr7+fr6fr10+fr7fr12  + fr8fr13,t4=fr1fr4+fr5fr8+fr6fr11+fr7fr13  + fr8fr14,t5=fr2fr2+fr6fr6+fr9fr9+fr10fr10  + fr11fr11,t6=fr2fr3+fr6fr7+fr9fr10+fr10fr12  + fr11fr13,t7=fr2fr4+fr6fr8+fr9fr11+fr10fr13  + fr11fr14,t8=fr3fr3+fr7fr7+fr10fr10+fr12fr12  + fr13fr13,t9=fr3fr4+fr7fr8+fr10fr11+fr12fr13  + fr13fr14,t10=fr4fr4+fr8fr8+fr11fr11+fr13fr13  + fr14fr14 , 
									 where frk denotes the filter response for monomial filter k. the first term relates to q1q1 , and the rest of the terms relate to q2q2 , in total q1q1 + q2q2 . if monomial filters are used instead of quadrature filters , the required number of 4d filters is thus decreased from 24 to 14 .
another advantage is that the monomial filters require a smaller spatial support , which makes it easier to preserve details and contrast in the processing .
a smaller spatial support also results in a lower number of filter coefficients , which decreases the processing time .
when the local structure tensor t has been estimated , it is then mapped to a control tensor c , by mapping the magnitude ( energy ) and the isotropy of the tensor .
the purpose of this mapping is to further improve the denoising . for 2d and 3d image denoising
, this mapping can be done by first calculating the eigenvalues and eigenvectors of the structure tensor in each element of the data .
the mapping is first described for 2d and then for 4d . in the 2d case ,
the magnitude 0 of the tensor is calculated as 
 ( 12)0=12+22 , 
									 where 1 and 2 are the two eigenvalues .
the magnitude 0 is normalized to vary between 0 and 1 and is then mapped to  with a so - called m - function according to 
 ( 13)=(0  0    +  +   ) , 
									 where  ,  , and  are parameters that are used to control the mapping .
the  variable is directly proportional to the signal - to - noise ( snr ) ratio of the data and acts as a soft noise threshold ,  mainly controls the overshoot ( that can be used for dynamic range compression or to amplify areas that have a magnitude slightly above the noise threshold ) , and  mainly controls the slope / softness of the curve . the purpose of this mapping is to control the general usage of highpass information .
the highpass information should only be used where there is a well - defined structure in the data .
if the magnitude of the structure tensor is low , one can assume that the neighbourhood only contains noise .
the isotropy 0 is in 2d calculated as 
 ( 14)0=21 
									 and is mapped to  with a so called mu - function according to 
 ( 15)=(0(1))(0(1))+((10)) , 
									 where  and  are parameters that are used to control the mapping ,  mainly controls the transition of the curve and  mainly controls the slope / softness .
the purpose of this mapping is to control the usage of highpass information in the nondominant direction , that is , the direction that is given by the eigenvector corresponding to the smallest eigenvalue .
this is done by making the tensor more isotropic if it is slightly isotropic , or making it even more anisotropic if it is anisotropic .
the control tensor c is finally calculated as 
 ( 16)c=e1e1t+e2e2 t , 
									 where e1 is the eigenvector corresponding to the largest eigenvalue 1 and e2 is the eigenvector corresponding to the smallest eigenvalue 2 .
the mapping thus preserves the eigensystem , but changes the eigenvalues and thereby the shape of the tensor . for matrices of size 2  2 and 3  3 , there are direct formulas for how to calculate the eigenvalues and eigenvectors , but for 4  4 matrices , there are no such formulas and this complicates the mapping .
it would of course be possible to calculate the eigenvalues and eigenvectors by other approaches , such as the power iteration algorithm , but this would be extremely time consuming as the mapping to the control tensor has to be done in each time voxel .
the mapping of the local structure tensor to the control tensor is in this work therefore performed in a way that does not explicitly need the calculation of eigenvalues and eigenvectors .
the tensor magnitude is first calculated as 
 ( 17)tmag=||t8||1/8 , 
									 where |||| denotes the frobenius norm .
the exponent will determine how close to 1 the estimated tensor magnitude will be ; a higher exponent will give better precision ,   but an exponent of 8 has proven to be sufficient in practice . to reduce the computational load ,
t is calculated as 
 ( 18)t2=tt , t4=t2t2,t8=t4t4 , 
									 where  denotes matrix multiplication .
0 is then calculated as 
 ( 19)0=tmagmax  ( tmag ) , 
									 where the max operator is for the entire data set , such that the maximum value of 0 will be 1 , 0 is then mapped to  by using the m - function . to map the isotropy ,
the structure tensor is first normalized as 
 ( 20)t^=ttmag , 
									 such that the tensor only carries information about the anisotropy ( shape ) .
the fact that t^ and i - t^ have the same eigensystem is used , such that the control tensor can be calculated as 
 ( 21)c=(i+(1)t^ ) , 
									 where i is the identity matrix .
the following formulas are an ad hoc modification of this basic idea , that do not explicitly need the calculation of the isotropy  and that give good results for our ct data .
the basic idea is that the ratio of the eigenvalues of the tensor change when the tensor is multiplied with itself a number of times , and thereby the shape of the tensor also changes .
this approach does not give exactly the same results as the original isotropy mapping , but it circumvents the explicit calculation of eigenvalues and eigenvectors .
a help variable t^f is first calculated as 
 ( 22)t^f = t^2(i+2(it^ ) ) , 
									 and then the control tensor c is calculated as 
 ( 23)c=  ( i(it^f)8(i+8t^f ) ) . 
  the resulting transfer function that maps each eigenvalue is given in figure 7 .
eigenvalues that are small become even smaller , and eigenvalues that are large become even larger . the result of this eigenvalue
eleven nonseparable reconstruction filters , one lowpass filter h0 of the zeroth order and 10 highpass filters h2,mn of the second order , with a spatial support of 11  11  11  11 time voxels are applied to the ct volumes .
the denoised 4d data i d is calculated as the sum of the lowpass - filtered data , ilp , and the highpass filtered data for each highpass - filter k , ihp(k ) , weighted with the components ck of the control tensor c , 
 ( 24)id = ilp+k=110ckihp(k ) . 
							
the result is that the 4d data is lowpass filtered in all directions and then highpass information is put back where there is a well - defined structure .
highpass information is put back in the dominant direction of the local neighbourhood ( given by the eigenvector related to the largest eigenvalue ) if the tensor magnitude is high .
highpass information is put back in the nondominant directions ( given by the eigenvectors related to the smaller eigenvalues ) if the tensor magnitude is high and the anisotropy is low .
all the processing steps of the denoising algorithm are given in table 2 . in our case
the ct data does not contain any significant structural information in the frequencies over /2 in the spatial dimensions , the volumes are therefore lowpass filtered and then downsampled a factor 2 in x , y , z. when the local structure tensor has been estimated , it is lowpass filtered , with a separable lowpass filter of size 5  5  5  3 , to improve the estimate in each time voxel and to make sure that the resulting reconstruction filter varies smoothly .
note that this smoothing does not decrease the resolution of the image data , but only the resolution of the tensor field .
after the tensor mapping , the control tensor is interpolated to the original resolution of the ct data . while the presented algorithm is straightforward to implement , spatial filtering with 11 reconstruction filters of size 11  11  11  11 ( 14  641 filter coefficients ) applied to a dataset of the resolution 512  512  445  20 requires about 375  000 billion multiplications .
this is the reason why the gpu is needed in order to do the 4d denoising in a reasonable amount of time .
one of the main drawbacks of the presented algorithm is that , using standard convolution , the number of valid elements in the z - direction ( i.e. , slices ) decreases rapidly .
if the algorithm is applied to a dataset of the resolution 512  512  34  20 , two slices are first lost due to the convolution with the lowpass filter of size 3  3  3 .
the monomial filters are of size 7  7  7  7 , thereby only 10 of the filter response slices are valid . during the lowpass filtering of each structure tensor component ,
another four slices are lost and then another four are lost during lowpass filtering of the control tensor .
the result is thus that only 2 valid slices are left after all the convolutions .
the same problem could exist in the time dimension , but since the heart cycle is periodic it is natural to use circular convolution in the time direction , and thereby all the time points are valid .
the loss of valid slices can be avoided by using normalized convolution , both for the lowpass filtering of the data before downsampling and the lowpass filtering of the tensor components . in normalized convolution , a certainty
a certainty - weighted filter response cwr is calculated as 
 ( 25)cwr=(cs)fcf , 
							 where c is the certainty , s is the signal , f is the filter ,  denotes pointwise multiplication , and  denotes convolution .
note that this simple version of normalized convolution ( normalized averaging ) can not be applied for the monomial filters and for the reconstruction filters , as these filters have both negative and positive coefficients .
it is possible to apply the full normalized convolution approach for these filters , but it will significantly increase the computational load .
  in this section , the gpu implementation of the denoising algorithm will be described .
the cuda ( compute unified device architecture ) programming language by nvidia , explained by kirk and hwu , has been used for the implementation .
the open computing language ( opencl )   could be a better choice , as it makes it possible to run the same code on any hardware .  
the cuda programming language can easily generate 2d indices for each thread , for example , by using algorithm 1 . to generate 3d indices
is harder , as each thread block can be three dimensional but the grid can only be two dimensional .
to generate 4d indices is even more difficult . to navigate in the 4d data ,
the 3d indexing approach described above is used , and the kernel is then called once for each time point . fast - fourier - transform ( fft- )
based filtering can be very efficient when large nonseparable filters of high dimension are to be applied to big datasets , but spatial filtering is generally faster if the filters are small or cartesian separable .
the main advantage with fft - based filtering is that the processing time is the same regardless of the spatial size of the filter .
the main disadvantage with fft - based filtering is however the memory requirements , as the filters need to be stored in the same resolution as the data , and also as a complex - valued number for each element . to see which kind of filtering that fits the gpu best , both spatial and fft - based filtering was therefore implemented . for filtering with the small separable lowpass filters ( which are applied before the data is downsampled and to smooth the tensor components ) ,
one easy way to implement 2d and 3d filtering on the gpu is to take advantage of the cache of the texture memory and put the filter kernel in constant memory .
the drawback with this approach is however that the implementation will be very limited by the memory bandwidth , and not by the computational performance .
another problem is that it is not possible to use 4d textures in the cuda programming language .
one would have to store the 4d data as one big 1d texture or as several 2d or 3d textures .
a better approach is to take advantage of the shared memory , which increased a factor 3 in size between the nvidia gtx 285 and the nvidia gtx 580 .
the data is first read into the shared memory and then the filter responses are calculated in parallel . by using the shared memory
, the threads can share the data in a very efficient way , which is beneficial as the filtering results for two neighbouring elements are calculated by mainly using the same data .
as multidimensional filters can be separable or nonseparable ( the monomial filters and the reconstruction filters are nonseparable , while the different lowpass filters are separable ) two different spatial filtering functions were implemented .
our separable 4d convolver is implemented by first doing the filtering for all the rows , then for all the columns , then for all the rods and finally for all the time points .
the data is first loaded into the shared memory and then the valid filter responses are calculated in parallel .
for the four kernels , 16  kb of shared memory is used such that 3 thread blocks can run in parallel on each multiprocessor on the nvidia gtx 580 .
the shared memory approach works rather well for nonseparable 2d filtering but not as well for nonseparable 3d and 4d filtering .
the size of the shared memory on the nvidia gtx 580 is 48  kb for each multiprocessor , and it is thereby only possible to , for example , fit 11  11  11  9 float values into it . if the 4d filter is of size 9  9  9  9 , only 3  3  3  1 = 27 valid filter responses can be generated for each multiprocessor .
a better approach for nonseparable filtering in 4d is to instead use an optimized 2d filtering kernel , and then accumulate the filter responses by summing over the other dimensions by calling the 2d filtering function for each slice and each time point of the filter .
our nonseparable 2d convolver first reads 64  64 pixels into the shared memory , then calculates the valid filter responses for all the 14 monomial filters or all the 11 reconstruction filters at the same time , and finally writes the results to global memory .
two versions of the convolver were implemented , one that maximally supports 7  7 filters and one that maximally supports 11  11 filters .
the first calculates 58  58 valid filter responses , and the second calculates 54  54 valid filter responses . as 64  64 float values only require 16  kb of memory , three thread blocks can run at the same time on each multiprocessor .
this results in 58  58  3 = 10092 and 54  54  3 = 8748 valid filter responses per multiprocessor . for optimal performance , the 2d filtering loop was completely unrolled by generating the code with a matlab script .
the 14 monomial filters are of size 7  7  7  7 , this would require 135 kb of memory to be stored as floats , but the constant memory is only 64  kb .
for this reason , 7  7 filter coefficients are stored at a time and are then updated for each time point and for each slice .
it would be possible to store 7  7  7 filter coefficients at a time , but by only storing 7  7 coefficients , the size of the filters ( 2.75  kb ) is small enough to always be in the cache of the constant memory ( 8  kb ) .
the same approach is used for the 11 reconstruction filters of size 11  11  11  11 .   while the cufft library by nvidia supports 1d , 2d , and 3d ffts , there is no direct support for 4d ffts .
as the fft is cartesian separable , it is however possible to do a 4d fft by applying four consecutive 1d ffts .
the cufft library supports launching a batch of 1d ffts , such that many 1d fft 's can run in parallel .
the batch of 1d ffts are applied along the first dimension in which the data is stored ( e.g. , along x if the data is stored as ( x , y , z , t ) ) . between each 1d fft , it is thereby necessary to change the order of the data ( e.g. , from ( x , y , z , t ) to ( y , z , t , x ) ) .
the drawback with this approach is that the time it takes to change order of the data can be longer than to actually perform the 1d fft .
the most recent version of the cufft library supports launching a batch of 2d fft 's . by applying two consecutive 2d fft 's
, it is sufficient to change the order of the data once , instead of three times .
a filter is padded with zeros to the same resolution as the data and is then transformed to the frequency domain .
to do the filtering , a complex - valued multiplication between the data and the filter is applied and then an inverse 4d fft is applied to the filter response . after the inverse transform ,   a fft shift is necessary ; there is however no such functionality in the cufft library .
when the tensor components and the denoised data are calculated , each of the four coordinates is shifted by using a help function , see algorithm 4 .  
as the monomial filters only have a real part or an imaginary part in the spatial domain , some additional time is saved by putting one monomial filter in the real part and another monomial filter in the imaginary part before the 4d fft is applied to the zero - padded filter .
when the complex multiplication is performed in the frequency domain , two filters are thus applied at the same time .
after the inverse 4d fft , the first filter response is extracted as the real part and second filter response is extracted as the imaginary part .
the main problem of implementing the 4d denoising algorithm on the gpu is the limited size of the global memory ( 3  gb in our case ) .
this is made even more difficult by the fact that the gpu driver can use as much as 100200  mb of the global memory .
storing all the ct data on the gpu at the same time is not possible , a single ct volume of the resolution 512  512  445 requires about 467  mb of memory if 32  bit floats are used .
storing the filter responses is even more problematic . to give an example , to store all the 11 reconstruction filter responses as floats for
a dataset of the size 512  512  445  20 would require about 103  gb of memory .
the denoising is therefore done for a number of slices ( e.g. , 16 or 32 ) at a time . for the spatial filtering ,
the algorithm is started with data of the resolution 512  512  51  20 and is downsampled to 256  256  26  20 .
the control tensor is calculated for 256  256  20  20 time voxels , and the denoised data is calculated for 512  512  39  20 time voxels . to process all the 445 slices
the algorithm is started with data of the resolution 512  512  31  20 and is downsampled to 256  256  16  20 .
the control tensor is then calculated for 256  256  10  20 time voxels , and the denoised data is calculated for 512  512  18  20 time voxels .
to process all the 445 slices requires 26 runs . to store the 10 components of the control tensor in the same resolution as the original ct data for one run with spatial filtering ( 512  512  39  20 ) would require about 12.2  gb of memory . as the control tensor needs to be interpolated a factor 2 in each spatial dimension , since it is estimated on downsampled data ,
another approach is used . interpolating the tensor is a perfect task for the gpu , due to the hardware support for linear interpolation .
the 10 tensor components , for one timepoint , are therefore stored in 10 textures and then the interpolation is done on the fly when the denoised data is calculated . by using this approach , only another 10 variables of the resolution 256  256  20 need to be stored at the same time . 

table 2 states the in and out resolution of the data , the used equations , and the memory consumption at each step of the denoising algorithm , for spatial filtering and fft - based filtering .
the out resolution refers to the resolution of the data that is valid after each processing step , as some data is regarded as non - valid after filtering operations .
the reason why the memory consumption is larger for the fft - based filtering is that the spatial filtering can be done for one slice or one volume at a time , while the fft - based filtering has to be applied to a sufficiently large number of slices and time points at the same time .
we were not able to use more than about 2  gb of memory for the fft - based filtering ; one reason for this might be that the cufft functions internally use temporary variables that use some of the memory .
since the source code for the cufft library is unavailable , it is hard to further investigate this hypothesis .
the 4d ct dataset that was used for testing our gpu implementation was collected with a siemens somatom definition flash dual - energy ct scanner at the center for medical image science and visualization ( cmiv ) .
the dataset contains almost 9000 dicom files and the resolution of the data is 512  512  445  20 time voxels . the spatial size of each voxel is 0.75  0.75  0.75  mm . during the image acquisition
the tube current is modulated over the cardiac cycle with reduced radiation exposure during the systolic heart phase . due to this , the amount of noise varies with time .
a comparison between the processing times for our gpu implementation and for a cpu implementation was made .
the used gpu was a nvidia gtx 580 , equipped with 512 processor cores and 3 gb of memory ( the nvidia gtx 580 is normally equipped with 1.5  gb of memory ) .
the used cpu was an intel xeon 2.4  ghz with 4 processor cores and 12  mb of l3 cache , 12  gb of memory was used .
, the openmp ( open multiprocessing ) library [ 38 , 39 ] was used , such that all the 4 processor cores work in parallel .
no other types of optimization for the cpu , such as sse2 , were used .
we are fully aware of the fact that it is possible to make a much better cpu implementation .
the purpose of this comparison is rather to give an indication of the performance of the cpu and the gpu .
if the cpu code would be vectorized , the cpu processing times can be divided by a factor 3 or 4 ( except for the fft which already is very optimized ) .
the processing times are given in tables 3 , 4 , 5 , and 6 .
the cpu processing times for the spatial filtering are estimates , since it takes several days to run the algorithm on the whole dataset .
the processing times for a multi - gpu implementation would scale rather linearly with the number of gpus , since each gpu can work on different subsets of slices in parallel . as our computer contains three gpus , all the processing times for the gpu can thereby be divided by a factor 3 .
to show the results of the 4d denoising , the original ct data was compared with the denoised data by applying volume rendering .
two volume renderers , one for the original data and one for the denoised data , run at the same time and were synced in terms of view angle and transfer function .  
figure 8 shows volume renderings of the original and the denoised data for different time points and view angels .
it is clear that a lot of noise is removed by the denoising , but since the denoising algorithm alters the histogram of the data , it is hard to make an objective comparison even if the same transfer function is applied . a movie where the original and the denoised data is explored with the two volume renderers was also made .
for this video , the data was downsampled a factor 2 in the spatial dimensions , in order to decrease the memory usage .
the video can be found at   http://www.youtube.com/watch?v=wflbt2sv34m .    by looking at the video ,
it is easy to see that the amount of noise in the original data varies with time .
the result is that 4d image denoising becomes practically possible if the gpu is used and thereby the clinical value increases significantly . to make a completely fair comparison between the cpu and
it has been debated   if the gpu speedups that have been reported in the literature are plausible or if they are the result of comparisons with unoptimized cpu implementations . in our opinion ,
the theoretical and practical processing performance that can be achieved for different hardware is not the only interesting topic . in a research environment ,
the ratio between the achievable processing performance and the time it takes to do the implementation is also important . from this perspective
, we think that our cpu - gpu comparison is rather fair , since about the same time was spent on doing the cpu and the gpu implementation . the cuda programming language was designed and developed for parallel calculations from the beginning , while different addons have been added to the c programming language to be able to do parallel calculations .
while it is rather easy to make the cpu implementation multithreaded , for example , by using the openmp library , more advanced cpu optimization is often more difficult to include and often requires assembler programming . while spatial filtering can be significantly slower than fft - based filtering for nonseparable filters , there are some advantages ( except for the lower memory usage ) .
one is that a region of interest ( roi ) can be selected for the denoising , compared to doing the denoising on the whole dataset .
another advantage is that filter networks [ 41 , 42 ] can be applied , such that the filter responses from many small filters are combined to the same filter response as from one large filter .
filter networks can reduce the number of multiplications as much as a factor 5 in 2d , 25 in 3d and 300 in 4d . to design and optimize a filter network
another problem is that the memory usage increases significantly when filter networks are used , since many filter responses need to be stored in memory .
filter networks on the gpu is a promising area for future research . from our results , it is clear that fft - based filtering is faster than spatial filtering for large nonseparable filters . for data sizes that are not a power of two in each dimension , the fft based
since medical doctors normally do not look at 3d or 4d data as volume renderings , but rather as 2d slices , the spatial filtering approach however has the advantage that the denoising can be done for a region of interest ( e.g. , a specific slice or volume ) .
it is a waste of time to enhance the parts of the data that are not used by the medical doctor .
the spatial filtering approach can also handle larger datasets than the fft - based approach , as it is sufficient to store the filter responses for one slice or one volume at a time .
recently , we acquired a ct data set with 100 time points , compared to 20 time points .
it is not possible to use the fft - based approach for this data set .
there are several reasons why the gpu speedup for the fft - based filtering is much smaller than the gpu speedup for the spatial filtering .
first , the cufft library does not include any direct support for 4d fft 's , and we had to implement our own 4d fft as two 2d fft 's that are applied after each other . between the 2d fft 's the storage order of the data
it can take a longer time to change the order of the data than to actually perform the fft .
if nvidia includes direct support for 4d fft 's in the cufft library , we are sure that their implementation would be much more efficient than ours .
second , the fft for the cpu is extremely optimized , as it is used in a lot of applications , and our convolver for the cpu is not fully optimized . the cuda programming language is only a few years old , and the gpu standard libraries are not as optimized as the cpu standard libraries .
[ 44 , 45 ] have created their own gpu fft library which has been proven to give better performance than the cufft library .
they circumvent the problem of changing the order of the data and thereby achieve an implementation that is much more efficient . in 2008 , their 3d fft was 5 - 6 times faster than the 3d fft in the cufft library .
third , due to the larger memory requirements of fft - based filtering it is not possible to achieve an as big speedup for the gpu implementation as for the cpu implementation .
if a gpu with a higher amount of global memory would have been used , the fft - based implementation would have been more efficient . as previously discussed in the paper ,
4d image processing in cuda is harder to implement than 2d and 3d image processing .
there are , for example , no 4d textures , no 4d ffts , and there is no direct support for 4d ( or 3d ) indices . however , since fmri data also is 4d , we have previously gained some experience on how to do 4d image processing with cuda [ 1820 ] .
the conclusions that we draw after implementing the 4d image denoising algorithm with the cuda programming language is thus that cuda is not perfectly suited for 4d image processing , but due to its flexibility , it was still possible to implement the algorithm rather easily
. it might seem impossible to have medical image data with more than 4 dimensions , but some work has been done on how to collect 5d data .
the five dimensions are the three spatial dimensions and two time dimensions , one for the breathing rhythm and one for the heart rhythm .
one major advantage with 5d data is that the patient can breathe normally during the data acquisition , while the patient has to hold its breath during collection of 4d data . with 5d data , it is possible to , for example , fixate the heart and only see the lungs moving , or fixate the lungs to only see the heart beating .
if the presented algorithm would be extended to 5d , it would be necessary to use a total of 20 monomial filters and 16 reconstruction filters . for a 5d dataset of the size 512  512  445  20  20 , the required number of multiplications for spatial filtering with the reconstruction filters would increase from 375 000 billion for 4d to about 119 million billion ( 1.19  10 ) for 5d .
the size of the reconstruction filter responses would increase from 103  gb for 4d to 2986  gb for 5d .
this is still only one dataset for one patient , and we expect that both the spatial and the temporal resolution of all medical imaging modalities will increase even further in the future . except for the 5 outer dimensions , it is also possible to collect data with more than one inner dimension .
this is ,   for example , the case if the blood flow of the heart is to be studied . for flow data ,
a three - dimensional vector needs to be stored in each time voxel , instead of a single intensity value .
to conclude , by using the gpu , true 4d image denoising becomes practically feasible .
our implementation can of course be applied to other modalities as well , such as ultrasound and mri , and not only to ct data .
the short processing time also makes it practically possible to further improve the denoising algorithm and to tune the parameters that are used .
the elapsed time between the development of practically feasible 2d   and 3d   image denoising techniques was about 10 years , from 3d to 4d the elapsed time was about 20 years . due to the rapid development of gpus ,
 the use of image denoising techniques is an important part of many medical imaging applications . 
 one common application is to improve the image quality of low - dose ( noisy ) computed tomography ( ct ) data .   
 while 3d image denoising previously has been applied to several volumes independently , there has not been much work done on true 4d image denoising , where the algorithm considers several volumes at the same time . the problem with 4d image denoising , 
 compared to 2d and 3d denoising , is that the computational complexity increases exponentially . in this paper 
 we describe a novel algorithm for true 4d image denoising , based on local adaptive filtering , and how to implement it on the graphics processing unit ( gpu ) . 
 the algorithm was applied to a 4d ct heart dataset of the resolution 512   512   445   20 . 
 the result is that the gpu can complete the denoising in about 25 minutes if spatial filtering is used and in about 8 minutes if fft - based filtering is used . 
 the cpu implementation requires several days of processing time for spatial filtering and about 50 minutes for fft - based filtering . 
 the short processing time increases the clinical value of true 4d image denoising significantly . 