volume visualization is an essential tool for exploring and analysing the anatomy of complex structures and phenomena .
it has been extensively used in various scientific and engineering arenas such as medical imaging , geoscience , microscopy , mechanical engineering , and others [ 13 ] .
the main two categories that have gained broad acceptance by scientific communities were volume and surface rendering .
each technique has its specific area of application that is associated with its advantages , but in contrast , it also has its disadvantages that give the opportunity for other techniques to survive [ 46 ] .
volume rendering has shown great significance on the visual interpretation of large amount of 3d scalar and vector data generated by multidimensional sensors , acquisition devices , and supercomputer simulations .
it concentrates on visualizing the desired internal features of volumetric objects and their bounding surfaces at the same time . in literature ,
several volume rendering algorithms have been presented either to improve the rendering speed of large datasets or to enhance the rendering quality of their reconstructed images . by and large ,
the rendering speed and quality are traded - off and no single complete algorithm that can deliver the optimum quality associated with maximum interactivity exists .
a radical domain - based categorization of volume rendering algorithms classifies them into spatial - domain and other - domain - based techniques such as frequency domain , compression domain or the wavelet domain .
domain - based methods operate by switching part of their computations to a different domain either to reduce the complexity of the rendering algorithm or to improve the performance of some operations that take considerable amount of time in the spatial domain . for a volume of size n , complexity of spatial - domain algorithms
is of order (n ) , since all the voxels composing this volume must be visited at least once to render a correct image .
although some algorithms use additional optimization techniques to reduce the number of traversed samples , this optimization is in general data dependent and subject to the size of the datasets .
this time - complexity limits the usage of spatial domain rendering algorithms for interactive environments in some applications . in such cases , frequency - domain - based techniques can be used alternatively .
frequency - domain volume rendering ( fdvr ) uses a 3d spectral representation of the volume to compute an image that looks like x - ray radiograph in (nlogn ) time relying on the projection - slice theorem .
then it reconstructs 2d projection at any viewing angle by resampling an extracted projection - slice along a perpendicular plane to the viewing direction followed by backtransforming this resampled slice to the spatial domain [ 8 , 9 ] .
obtaining the spectral representation of the volume is the most computationally intensive step in the algorithm due to its (nlogn ) complexity .
nevertheless , this rendering algorithm is extremely efficient because this step is executed once in a preprocessing stage .
fdvr is a generic technique that can work with any frequency transform to switch between the spatial and frequency domains . our rendering pipeline will be based on fourier transform and so , the rendering method can be called fourier volume rendering .
section 4 elaborates the results gained by mapping the entire rendering pipeline to run on the gpu and section 5 concludes the paper .
although the technique was introduced around 20 years ago , the literature behind fdvr in general and fvr in particular is quite scarce . in this section
, we will try to summarize the most notable contributions that have been presented since that time in a nutshell . in 1992 , dunne et al .
malzbender systematically extended the projection - slice theorem to 3d to be a basis for volume rendering .
based on fft , he presented an implementation of the basic fvr pipeline and briefly discussed some considerations of resampling the frequency domain and their significant effects on the reconstruction quality of the generated projection image .
he proposed to use a hamming - windowed sinc reconstruction filter to reduce the aliasing accompanied with resampling the frequency spectrum .
for the same mission , grosso and ertl proposed alternatively biorthognal wavelet reconstruction filters . as a drawback ,
the basic fvr integral does not provide a direct way of modelling emission and scattering in the participating medium .
this limitation was a fundamental disadvantage associated with frequency domain rendering that is reflected as lack of occlusion in the resulting reconstructions .
levoy has partially restored the lost visual cues in the basic pipeline by applying several shading models that are linear combinations of fourier projections .
his extension included depth cueing , directional shading , and lambertian reflection with spherical illumination .
the result of levoy 's work did not exhibit real occlusion , yet it provided acceptable depth and shape cues that can simulate the existence of the missing occlusion .
one significant disadvantage that limited the usage of his implementation was the demand of several copies of the volume , which consequently imposes high memory requirements for initial preprocessing . in cooperation with totsuka ,
the same results have been obtained after considering alternative frequency domain methods by processing the frequency response of the volume data , which dramatically reduced the memory required before . until that time
, fvr did not gain that broad acceptance by the scientific visualization and medical communities due to the lack of illumination models in the fourier domain . in 2002 , entezari et al . enhanced the projections quality by incorporating various illumination models into the fvr pipeline .
the first model was based on gamma - corrected hemispherical shading that was proposed by scoggins et al .   and
fvr suffered from reduced reconstruction quality that limited its usefulness in particular medical applications . in , a solution for enhancing fvr using contour extraction was proposed .
it provided a flexible method for extracting material boundaries on surfaces in the fourier space .
additionally , it included enhancement of several features for revealing important spatial relationships between interior and exterior structures making it an attractive tool for improved x - ray - like investigations for a given dataset .
cheng and ching   also presented various methods for designing fvr transfer function based on bezier curves and b - splines .
partially accelerated the rendering pipeline on the gpu by mapping the split - stream - fft on the gpu .
as a major tool in digital radiography , ntasis et al .   provided a web - based fourier volume renderer for real - time preview of digital reconstructed radiographs ( drrs ) .
their implementation examined carefully the resampling issues of the frequency domain to generate remote high quality drrs . extending the technique beyond operating on volume data with regular grids ,
this framework abstracts the complex implementation details of the rendering pipeline to allow imaging researchers to develop image enhancement techniques without having prior knowledge of the opengl pipeline .
have implemented a fvr heterogeneous rendering pipeline that uses the gpu to execute the rendering stage relying on fragment shaders . however , this approach is limited if the current capabilities of unified computing gpu architectures are considered
. practical implementation of fvr is complicated by two main factors which arise when the projection - slice theorem is applied to discrete sampled data .
first of all , conventional fft algorithms yield frequency - domain output data which is not ideally structured for resampling . to mitigate this effect ,
it is necessary to add high performance multidimensional fft - shift stages to the rendering pipeline to rearrange the data .
moreover , sampling the frequency domain is equivalent to the replication of the signal in the spatial domain that ultimately accounts for the appearance of ghosting artifacts in the reconstructed images .
this issue is resolved by zero - padding the volume in the spatial domain and using high - order interpolation filters in the frequency domain .
the plain fvr algorithm is briefly illustrated to simplify the explanation of the contextual classification of the rendering pipeline . based on this classification ,
an efficient strategy has been considered to map this pipeline to run entirely on the gpu .
the spatial volume must be shifted by a 3d fft - shift operation to set its center at the origin of the 3d space .
the frequency spectrum of this volume is then obtained by a forward 3d fft operation .
consequently , another 3d fft - shift operation is required to center its zero - frequency component to prepare it for correct slice extraction .
afterwards , a 2d projection - slice passing through the origin of the spectral volume , with a normal that is parallel to the viewing direction , is carefully extracted .
this slice represents the 2d fft of the desired projection and thus , the reconstructed image can be directly obtained by a 2d inverse fft operation . to reduce the aliasing and ghosting artifacts in the reconstructed image ,
after resampling the extracted slice , it is backtransformed to the spatial domain by an inverse 2d fft operation . the resulting image from this inverse transformation
is shifted by another 2d fft - shift operation to move the center of the image from the edge to origin of the grid used to display the image .
changing the viewing angle implies rotating the 3d spectrum to extract a new projection - slice that corresponds to this angle .
the algorithm sequence is graphically depicted in figure 1 . according to the sequence of operations in this algorithm ,
the pipeline could be divided into two consecutive stages : a preprocessing stage and a rendering loop .
the preprocessing stage is executed only once to prepare the 3d spectrum of the input spatial volume .
the rendering loop is running continuously to generate different projections according to the input viewing angle .
the main function of the preprocessing stage is limited to loading a volume of interest , obtaining its frequency spectrum and preparing it for slicing .
after the preparation of the spectral volume , the rendering loop is executed to generate different projection images by extracting a projection - slice according to the input viewing angle , then resampling it , and finally backtransforming it to the spatial domain to generate the reconstructed image . from another perspective , the pipeline could be split according to functionality into two complementary cores or contexts : a computational context and a rendering context .
cooperatively , both contexts complement each other functional - wise , but from the sequence point of view , they are overlapping and one must switch from one context to the other to flow through the pipeline . building the entire fvr pipeline on the gpu
in general , an adequate way of constructing gpu - based pipelines is to have a nave and valid reference implementation of this pipeline on the cpu in advance .
having this cpu - based pipeline tested and validated , every stage of it can be afterwards implemented independently on the gpu and then validated to ensure similar results from the cpu - based implementation until getting the pipeline working integrally on the gpu .
this implementation strategy is a rule of thumb to increases flexibility and efficiency for designing and building gpu pipelines .
according to the aforementioned contextual functional classification of the fvr pipeline , the cpu is employed for executing the computational context while the gpu can be used to implement the rendering contexts relying on any graphics application programming interface ( api ) like opengl or direct 3d .
this heterogeneous implementation is easier to handle than a full cpu - based pipeline because the 3d spectral volume can be efficiently represented by 3d opengl textures .
this efficiency comes from their hardware - acceleration support and built - in interpolation schemes compared to using a 3d array on the cpu and doing the interpolation step manually .
the reference hybrid implementation of the fvr pipeline has been adopted from a previous single threaded one discussed in .
this volume is then spatially shifted and its frequency spectrum is obtained by a forward 3d fft operation . after activating the first opengl off - screen rendering context , this spectral volume is mapped to a 3d texture that is allocated on and uploaded to the gpu memory .
this texture is intersected with a proxy quadrant passing through the origin of the frequency domain . the resulting projection - slice
is then directed to a frame buffer object ( fbo ) and packed in a 2-component 2d texture attached to this fbo .
the contexts are then switched to move back to the computational context on the cpu , where the extracted projection - slice is downloaded from the 2d texture .
the communication between the computational and off - screen rendering contexts is shown in figure 3 .
afterwards , the computational context is resumed by resampling the projection - slice to remove the ghosting artifacts .
the resampled slice is then back - transformed by an inverse 2d fft operation to the spatial domain to produce a shifted image of the projection .
the final reconstructed image is then packed into a 2d opengl texture and uploaded via the command gltex2d to reside on the gpu memory .
finally , the on - screen rendering context is activated and the texture that contains the final projection image is sent to the frame buffer to be displayed .
as usual , each nave approach has its accompanied bottlenecks that have to be investigated for either their complete removal or at least minimizing their performance overheads .
the hybrid implementation of this technique obviously lacked the soul of interactivity due to several bottlenecks . in this section
, the hybrid pipeline is analyzed to suppress its main bottlenecks and also to optimize its flow in order to get ready to have it entirely mapped to the gpu with both of its computational and rendering contexts . a main concern that significantly affects the performance in the computational context is the fft operations .
the 3d fft operation takes a considerable amount of time if the size of the input volume was 64 or more .
although it is executed once during the pre - processing stage , performing this operation on several volumes to reveal any depth cues will introduce a real bottleneck and this will consequently elongate the pre - processing stage . additionally , the inverse 2d fft operation is executed on a per - frame basis . in turn , reducing the time consumed by this operation will be reflected as an order of magnitude enhancement in the frame rate .
the performance of the rendering loop is affected by the resampling stage , which represents the most critical bottleneck in this pipeline .
two of them are used to iterate on each dimension of the slice and the other two are used for the filter kernel . for a projection - slice of size 512 , the resampling operation takes a considerable amount of time that eventually blocks a real - time rendering loop . additionally , executing the computational context on the cpu and the rendering ones on the gpu
this transfer limits the overall performance of the pipeline to the bandwidth of the communication channel between the cpu and the gpu , which is much less than that on the gpu internally . mapping the fvr pipeline
to be entirely running on the gpu is a significant step that will leverage the overall performance of the reconstruction process .
although programming the gpu for executing nongraphics algorithms has been simplified by the introduction of generic gpu apis like cuda and opencl , but it is far from trivial to effectively exploit this titanic computing power provided by the gpu .
this point makes the communication between opengl contexts and the cpu comparatively simple . in the gpu - based approach
, the implementation of opengl contexts remains the same as that of the previous one , while the computational core of the pipeline will be implemented in a cuda context .
however , both contexts run on the gpu , but the interoperability between them takes a different mechanism and some workarounds for successful and efficient communication . the issues behind opengl interoperability with cuda are discussed with some sort of details in [ 31 , 32 ] . the first stage for mapping the computational core of the pipeline into a cuda context is writing analogous device kernels to the respective c functions of the cpu - based computational context . in our mapping , the following kernels have been written.fft_shift_3d_real , which wraps around the 3d spatial volume.fft_shift_3d_complex , which wraps around the 3d complex spectrum.fft_shift_2d_real , which wraps around the resulting image from the inverse 2d fft operation.resample_slice , which executes the high - order resampling operation.repack_array , which replaces the complex array resulting from the 3d cufft operation by another 1d alternative to match the format of the opengl 3d spectrum texture .
repack_array , which replaces the complex array resulting from the 3d cufft operation by another 1d alternative to match the format of the opengl 3d spectrum texture .
the fft - shift kernels have been adopted from the proposed implementations in [ 33 , 34 ] .
the following fft kernels were designed to encapsulate the fft plans that were originally implemented within the cufft library .cuda_fft_3d , which executes the forward 3d fft operation.cuda_fft_2d , which executes the inverse 2d fft operation .
cuda_fft_2d , which executes the inverse 2d fft operation . in this gpu - accelerated pipeline , shown in figure 5 ,
the cpu is only used for loading the input volume and controlling the flow of the pipeline .
this control includes switching between the different contexts , copying data from the cpu to the gpu , and dispatching opengl commands .
nevertheless , the entire flow of the pipeline is started and terminated on the gpu .
the opengl rendering contexts employed in the hybrid pipeline will be reused in this implementation to interoperate with the cuda computational context .
however , there will be slight modifications concerned with the communication mechanisms between the different contexts .
once the volume dataset of interest is loaded to the cpu side , it is directly shoved to reside on the gpu global memory by a cudamemcpy operation . on the gpu side ,
cuda_fft_3d function is then invoked to obtain the 3d frequency spectrum of the input spatial volume .
afterwards , the fft_shift_3d_complex kernel is invoked to set the center of the spectral volume at the origin of the k - space .
this array will be mapped to the opengl off - screen rendering context for performing the projection - slice extraction operation .
this complex array is converted to an alternative one that is compatible with an opengl 3d texture with two components using the repack_array kernel .
although the origin - aligned spectral volume array is located on the gpu memory , it can not be accessed directly by opengl because it is located in the cuda address space .
this array can be read back to the cpu and then uploaded again to the gpu to reside on the texture memory of the cuda context .
this memory transfer operation is extremely useless and can be solved by using an opengl pixel buffer object ( pbo ) that is registered in the cuda memory space .
this workaround makes it possible to directly copy the spectral array to a buffer that is shared between the different contexts of the pipeline .
after mapping the spectral texture to opengl memory space , we have to switch contexts to load the opengl off - screen rendering context .
the 3d spectral texture is bound and intersected with a polygon , with the same sizes of the reconstructed image , and the result of this intersection is directed to the fbo and stored in the attached 2d texture . by projecting the desired slice in this fbo , the opengl off - screen rendering context is terminated and the rendering pipeline switches back to the cuda context .
a reference to the attached texture to the fbo is mapped to the cuda context using cuda array in order to be accessible for subsequent kernel calls .
this cuda array that carries the projection - slice is then mapped to a cuda texture object .
the resample_slice kernel is then invoked , so the resampling operation is performed on this texture object and the results are stored directly into another 1d array with compatible format with the next cufft operation .
afterwards , the cuda_fft_2d function is executed on the resampled slice to generate the projection image .
this resulting image is shifted , and thus , an fft - shift operation is considered by invoking the fft_shift_2d_real kernel . this correct image is ready to be displayed , but it has to be mapped in advance to an opengl 2d texture .
this requires mapping the resulting array from the previous fft - shift operation to an opengl texture using a pbo . after this mapping operation
, the on - scree opengl rendering context is activated and the final image is pushed to the frame buffer for display .
changing the viewing angle repeats this flow by switching back to the opengl off - screen context to extract another projection - slice going through all the subsequent stages over and over until the termination of the running process .
although the opengl off - screen context is implemented the same way as it was done in the hybrid pipeline , the extracted projection - slice in the frame buffer object can not be read back by the cpu any more .
alternatively , it is mapped directly to the cuda memory space by creating another pbo and registering this buffer object with the cuda memory space .
this allows direct attachment of the extracted slice data from the opengl off - screen context to be accessible and callable from the different cuda kernels . due to the optimization of the texture caches for 2d spatial locality , reading device memory via texture fetches
is more performing than reading from the global memory . in that essence and to exploit the low - latency memory access associated with the texture memory , the extracted projection - slice is directly packed in a cuda array .
this array can be easily mapped to a cuda texture which is writable from the opengl context via the cuda array .
this texture leads to a significant speed - up in the resampling stage compared to the nave for loop implementation in the hybrid pipeline .
the context switching operations were accelerated due to the direct connection between the different cuda and opengl contexts relying on the shared buffer objects .
these buffer objects are registered within all the contexts to allow direct mapping of the resulting textures from the computational context to the rendering ones .
this permits efficient data sharing between the different contexts without any memory copies between the cpu and the cpu at all .

 figure 8 shows the resulting radiographs of four medical datasets with different sizes and three intensity scaling factors .
all the datasets are organized in regular 3d cartesian grids [ 37 , 38 ] . to show the artifacts associated with resampling the extracted projection slice ,
three different interpolation schemes have been applied : point - sampling , trilinear interpolation , and windowed - sinc interpolation .  
figure 9 reflects how the selection of the interpolation scheme can significantly affect the rendering quality of the projection image .
the projection reconstructed in figure 9(a ) does not exhibit any ghosts because it is orthogonal and all the replicas are not present in the view . in figure 9(b ) , the ghosting artifacts are very apparent due to the usage of the point sampling filter , which clamps the missing sample to the nearest available value . in figure 9(c ) , the intensities of the replicas are reduced when trilinear interpolation scheme is applied . in fact , the nearest - neighbour and tri - linear interpolation schemes have almost the same performance for resampling projection - slices with sizes less than or equal to 512 , but for larger volumes , the performance cost of trilinear interpolation is a bit high . however , this overhead is not significantly degrading the performance of the entire pipeline . 

figure 10 shows the result of zero - padding the spatial volume to suppress the overlapping between the central projection image and the surrounding replicas in the spatial domain .
it has to be noticed that the zero - padding operation comes with no overhead on the performance although it dramatically increases the memory requirements by a factor of 8 for 100% zero - padding . to minimize the artifacts associated with the nature of the technique
, the rendering pipeline should afford a high order interpolation filter combined with a preprocessing zero - padding operation .
a practical filter that can perform this operation is a hamming windowed - sinc interpolation filter with width 5 .  
figure 11 shows the result of combining this reconstruction filter with zero - padding the spatial volume on an oblique projection of the visible male datasets . on the performance side and in order to highlight the improvements gained by porting the computational context of the pipeline from the cpu to an alternative cuda context , we have analysed the profiling results for every stage individually and then we demonstrate and accounted for the overall speed - ups gained for the entire pipeline .
the rendering pipeline has been benchmarked on a workstation shipped with an intel core i7 - 4770 cpu running at 3.4  ghz with 8  mbyte of cache and 12  gbytes of ddr3 memory .
the application was compiled on ubuntu 14.04 using gcc 4.7.3 and nvidia cuda compiler nvcc 6.0.1 .
the first one was an nvidia geforce gt 640 , which is considered a midrange commodity gpu that costs now almost $ 100 .
this one is much more powerful than the first gpu and can be ranked as a state - of - the - art one since it has 1536  cuda cores and four gbytes of gddr5 memory .
it will be referred later in the text and the benchmarks by gpu 2 . additionally , and on the software side ,
the benchmarks were generated for different volumes of sizes 128 , 256 , and 512 .
although higher speed - ups can be achieved if the volume size goes beyond this limit , unfortunately , larger volumes can not fit in the memory of any of the employed gpus .  
figure 12 shows the benchmarking results for every stage in the pipeline to elaborate the difference between the cpu performance in comparison to the two gpus . 

table 1 summarizes the average execution times for the different stages of the pipeline and aggregates the entire pipeline performance for a volume dataset of size 512 .
we also reference the hybrid implementation to relatively show that our implementation has gained a speed - up of 117x for the entire pipeline and 29x for the rendering loop .
the gpu - accelerated pipeline requires uploading the spatial volume to the device global memory .
this step takes on average less than 100 and 40 milli - seconds for a volume of size 512 on gpu 1 and gpu 2 , respectively .
thanks to the interoperability mechanisms between cuda and opengl , there are no further data transfer operations between the different stages of the rendering pipeline . these profiles can be reduced if compared to a multithreaded implementation . however , building a multithreaded opengl application is quite complex , and thus a pure gpu implementation is much more appropriate to consider rather than constructing a multithreaded application for benchmarks comparison .
this paper presented a high performance implementation of the fourier volume rendering algorithm on cuda - capable gpus .
the literature of fvr was briefly covered and followed by a detailed explanation of the plain rendering algorithm by introducing the notions of computational and rendering contexts .
our implementation strategy considered in advance building a reference hybrid rendering solution where the computational context is executed on the cpu and the rendering one on the gpu using opengl .
this reference pipeline was then imported in a step - wise fashion to run entirely on the gpu by mapping the computational context to run within a cuda context and redesigning the rendering contexts to communicate between cuda and opengl to optimize the performance of the rendering loop . using a 512 dataset ,
the pure gpu implementation outperformed the hybrid one by a factor of 117x for the entire pipeline and nearly 29x of speed up for the rendering loop .
the gpu - accelerated pipeline will be extended in the future to include depth cues and different shading models that can improve the visual appearance of the resulting digital radiographs . a client - server distributed model of the pipeline
will be also considered to allow rendering large - scale datasets on multiple computing nodes .
 fourier volume rendering ( fvr ) is a significant visualization technique that has been used widely in digital radiography . as a result of its (n2logn ) time complexity 
 , it provides a faster alternative to spatial domain volume rendering algorithms that are (n3 ) computationally complex . relying on the fourier projection - slice theorem , this technique operates on the spectral representation of a 3d volume instead of processing its spatial representation to generate attenuation - only projections that look like x - ray radiographs . due to the rapid evolution of its underlying architecture , 
 the graphics processing unit ( gpu ) became an attractive competent platform that can deliver giant computational raw power compared to the central processing unit ( cpu ) on a per - dollar - basis . 
 the introduction of the compute unified device architecture ( cuda ) technology enables embarrassingly - parallel algorithms to run efficiently on cuda - capable gpu architectures . in this work , 
 a high performance gpu - accelerated implementation of the fvr pipeline on cuda - enabled gpus is presented . 
 this proposed implementation can achieve a speed - up of 117x compared to a single - threaded hybrid implementation that uses the cpu and gpu together by taking advantage of executing the rendering pipeline entirely on recent gpu architectures . 
1. Introduction
2. Fourier Volume Rendering Literature
3. Algorithm and Implementation
4. Results and Discussion
5. Conclusion and Future Work