/articles/toma

 X-ray tomography is a powerful tool to investigate materials and small animals at micro- and nano-scale~\cite{withers2007nanotomography}. Information about X-ray attenuation or/and phase changes in the sample are used to reconstruct its internal structure. At synchrotron imaging beamlines, samples are typically placed at a rotating stage in front of  a pixel detector. It registers a series of 2D projections while the sample is turning. Various imaging techniques are then used to reconstruct the full 3D volume from the series of projections. Recent advances in X-ray optics and detector technology have paved the way for a variety of new X-ray imaging experiments aiming to study dynamic processes in materials and to analyze small organisms in vivo. At the Swiss Light Source (SLS) scientists were able to take high quality 3D snapshots of 150 Hz oscillations of a blowfly flight motor~\cite{mokso2015fly}. A temporal resolution of 20 milliseconds was achieved during a stencil test performed at SLS~\cite{maire2016} and also in the analysis of morphological dynamics of fast-moving weevils at the ANKA synchrotron at KIT~\cite{dosSantosRolo2013invivo}.

 

 To achieve these results, the instrumentation used at imaging beamlines has recently undergone a major update. The nowadays installed streaming cameras are able to deliver up to hundreds of thousands of frames per second with a continuous data rate up to 8 GB/s~\cite{marone2017}. Newly developed control systems at ANKA~\cite{vogelgesang2016}, SLS~\cite{marone2017}, and other synchrotron facilities use the acquired imaging information to track the processes under study and adjust the instrumentation accordingly. These control systems rely highly on the performance of the integrated image processing framework. 

 

 Faster acquisition and a  higher level of automation is essential to study dynamic phenomena, but also enables experiments with significantly increased the sample throughput. For example, in 2015 Diamond Light Source (DLS) has reported that typically about 3000 scans are recorded during 5 days of operation at a single imaging beamline~\cite{atwood2015}. Consequently, the amount of data generated at imaging beamlines quickly grows and results in a steep rise in the computational demand of these techniques. In order to achieve higher temporal resolution and prolong experiment duration, advanced methods are developed that incorporate a priori knowledge in the reconstruction procedure. These methods are able to produce high quality images from undersampled and underexposed measurements, as demonstrated by~\cite{mirone2014dict, eyndhoven2015}. Unfortunately these methods are computationally significantly more demanding then traditional reconstruction algorithms and further increase the load to the computing infrastructure~\cite{ashkarin2015}.

 To tackle the performance challenge several reconstruction frameworks have been developed and optimized to use parallel capabilities of nowadays computing architectures. At SLS \emph{GridRec}, a fast reconstruction approach optimized for conventional CPU technology, has been adopted~\cite{marone2012regridding}. The reconstruction is scheduled across a dedicated cluster and reconstructs a 3D image within a couple of minutes~\cite{marone2017}. Other frameworks use GPUs to accelerate the computation and are able to achieve minute-scale reconstructions at a single node equipped with multiple GPU adapters. PyHST is developed at ESRF and uses the CUDA framework to offload image reconstruction to NVIDIA GPUs~\cite{chilingaryan2011gpu}. The second version of PyHST provides also a number of iterative reconstruction techniques~\cite{mirone2014pyhst2}. The UFO parallel computing framework is used at ANKA synchrotron to realize in-vivo tomography and laminography experiments~\cite{vogelgesang2012ufo,vogelgesang2016dgma}. It constructs a data processing workflow by combining basic building blocks in a graph structure. OpenCL is used to execute the reconstruction at parallel accelerators with a primary focus on NVIDIA and AMD GPUs. ASTRA is a fast and flexible development platform for tomographic algorithms with MATLAB and python interfaces~\cite{aarle2016astra,palenstijn2017astra}. It is implemented in C++ and uses CUDA to offload computations to GPU. Several other frameworks are based on ASTRA libraries to provide GPU-accelerated reconstruction, for instance the Savu framework at DLS~\cite{atwood2015} or TomoPy at the Advanced Photon Source (APS)~\cite{gursoy2014tomopy}. Recent versions of TomoPy also support \emph{UFO} and \emph{GridRec} as backends.  All of the GPU-accelerated frameworks are capable to distribute the computation across the cluster of computation nodes as well.

 

 While most of the nowadays imaging frameworks are heavily rely on the parallel hardware to speed-up reconstruction, the implementations are rarely optimized for the specific architectures but utilizing the most common denominator instead. On other hand, the hardware architectures differ significantly~\cite{zhang2014performance}. Organization of memory and cache hierarchies, performance balance between different types of operations, and even type of parallelism varies. A significant speed-up is possible if the details of the specific architecture are taken into the consideration as illustrated in~\cite{csa2011sc}. This is especially important if the reconstruction is running as part of control workflow. A reduced latency may allow to track faster processes and to improve the achieved spatial and temporal resolution. Due to imminent communication latencies, the same result can not always be reached by scaling the reconstruction cluster.

 

 For online monitoring and control, normally fast analytical methods are used to reconstruct 3D images. There are two main approaches, Filtered Back Projection (FBP) and methods based on the Fourier Slice Theorem~\cite{natterer2001mathematical}. The later methods are asymptotically faster, but due to the involved interpolation in the Fourier domain are more sensitive to the quality of the available projections. For typical geometries Fourier-based methods are several times faster using the same computing hardware~\cite{rshkarin2015} and should be preferred if the computing infrastructure is limited to general-purpose processors only~\cite{marone2017}. Recent study suggests to implement back projection as convolution in log-polar coordinates in order to gain high reconstruction speed with interpolation in the image domain~\cite{anderson2016fbpconv}. However, this new method has not yet been adopted in production environments. Still, Filtered Back Projection is the method of choice, largely due to it simplicity and robustness. Therefore, the efficiency of the FBP implementation is still crucial for the currently used real-time monitoring and control systems. Furthermore, since methods used for low dose tomography normally consist of sequences of forward and back projections, a faster implementation of the back projection  lower the computational demands for high-quality offline reconstruction and might reduce the required hardware investments.

 In this article we will review differences between recent parallel architectures and explain how the FBP algorithm can be mapped on the existing parallel hardware efficiently. We will focus on the architecture-specific optimizations and will only briefly mention the organization of data flow as it is already explained in the literature~\cite{chilingaryan2011gpu,vogelgesang2016dgma}. We also do not cover scaling issues since the proposed optimizations can be easily integrated in existing frameworks like ASTRA, PyHST, or UFO which are already providing multi-GPU and GPU-cluster support. The article is organized as follows. The hardware setup, software configuration, and pseudo-code conventions are listed in \sectionname~\ref{section:setup}. The short introduction to parallel architectures required to understand the proposed optimizations is given in \sectionname~\ref{section:arch}. In this section we also highlight the differences between considered parallel architectures. Filtered Back Projection algorithm and the state-of-the-art implementation are presented in \sectionname~\ref{section:tomo}. A number of optimizations to the state-of-the-art implementation of the back-projection algorithm are proposed in \sectionname~\ref{section:texrec}. An alternative implementation relaying on a different set of hardware resources is developed in \sectionname~\ref{section:alurec}. A hybrid approach combining both approaches to fully utilize all hardware resources is presented in \sectionname~\ref{section:hybridrec}.  The achieved performance improvements are discussed in \sectionname~\ref{section:summary}.