Multi-sheet surface rebinning methods for reconstruction from asymmetrically truncated cone beam projections: II. Axial deconvolution

In airport baggage scanning it is desirable to have a system that can scan baggage moving at standard conveyor belt speeds. One way to achieve this is to use multiple electronically switched sources rather than a single source on a mechanically rotated gantry. In such a system placing the detectors opposite the sources would obstruct the beam, so they have to be offset (hence offset multi-source geometry). This results in asymmetrical axial truncation of the cone beam projections. As such projections do not constitute complete data in the sense of integral geometry, the standard cone beam reconstruction algorithms do not apply. In this series of papers we introduce a new family of rebinning methods for reconstruction from axially asymmetrically truncated cone beam projections. In the first paper we discussed the approximation of the data on the multi-sheet surface with the truncated projection data obtained from offset multi-source geometries. In this second paper we focus on the recovery of the volumetric image from the reconstruction of data rebinned to multi-sheet surfaces. Multi-sheet rebinning effects an implicit relation between the fan beam transforms on the individual sheets and the rebinned data. This relation in conjunction with the linearity of the ray transform allows us to formulate the deconvolution problem for the recovery of the volume from a stack of reconstructed images on multi-sheet surfaces. We discuss the errors in the right-hand side of the deconvolution problem (reconstruction on multi-sheet surfaces) resulting from rebinning approximation. We introduce convolution matrix models based on the distribution of the distances of the rays from the multi-sheet surface, which considerably improve the data model fit and in turn lead to a superior solution. Multiple strategies for solution of the deconvolution problem are discussed and an efficient and robust implementation is presented, which makes the method capable of real time reconstruction. We conclude with some reconstruction results from simulated as well as real data collected with a Rapiscan RTT80 scanner.

some reconstruction results from simulated as well as real data collected with a Rapiscan RTT80 scanner.
(Some figures may appear in colour only in the online journal)

Introduction
The demand for fast (near real time) fully three-dimensional x-ray imagery motivated the development of a new scanning geometry offset multi-source geometry. In standard cone beam scanners, the detector is placed centrally opposite the source of radiation, and the sourcedetector assembly rotates around the object while the object is translated resulting in a helical trajectory of the source relative to the object. In contrast in the offset multi-source geometry the mechanical motion of the gantry was eliminated by introduction of a static ring of electronically switchable sources and a static cylindrical array of detectors, which is axially offset from the ring of sources, so that the (back of the) detector array does not obstruct the x-ray cone emitted from the active source. As a consequence the measured cone beam projections are axially asymmetrically truncated i.e. data is collected only in a subset of the Tam-Danielsson window, rendering standard cone beam reconstruction techniques not applicable.
In this series of papers we propose a new family of multi-sheet rebinning methods, for reconstruction from the axially truncated projections. Multi-sheet rebinning methods [1,3] combine analytical and numerical techniques to provide efficient methods with local data dependence, as is the case for single-sheet surface rebinning methods [4-7, 9-12, 16, 19, 20], while allowing for high resolution reconstruction from, in the sense of integral geometry, severely limited data.
Multi-sheet rebinning methods comprise three stages. First the truncated cone beam projections are rebinned to a stack of two-dimensional (2D) problems, each corresponding to fan beam transform on a multi-sheet surface. In the second stage the fan beam data on multisheet surfaces is reconstructed using full scan filtered backprojection algorithm. Finally, the reconstructed data on a stack of multi-sheet surfaces is deconvolved to provide the volumetric image, as opposed to simple interpolation used in single-sheet surface rebinning methods. The approximation of the data on multi-sheet surfaces and the problem of finding an optimal multi-sheet surface and rebinning function pair was the subject of the first paper in this series [3]. This second paper deals with the axial deconvolution problem.
The remainder of the paper is organized as follows. In the next section we briefly summarize the results on rebinning to multi-sheet surfaces obtained in [3]. In section 3 we discuss reconstruction of the fan beam data on a multi-sheet surface. The impact of the rebinning errors is discussed in section 4. In section 5 we set up the deconvolution problem for the recovery of the volumetric image from a stack of reconstructed images on the multi-sheet surfaces. We introduce a realistic convolution matrix model based on the back projection formulae and demonstrate its effectiveness compared to a simple model assuming that the rebinning of the rays to the multi-sheet surface is exact. In section 6 we discuss an efficient implementation of the deconvolution step. Different solution strategies for the deconvolution problem are discussed in section 7. Section 8 demonstrates the performance of the multi-sheet surface rebinning methods on simulated phantom data as well as the real data acquired with a Rapiscan Systems RTT80 scanner (RTT). We conclude with a summary of the multi-sheet surface rebinning methods and some prospective research directions.

Multi-sheet surface rebinning
Offset multi-source geometry as deployed in the Rapiscan RTT80 scanner features a static ring of sources, and multiple static rings of detectors located in axially offset planes parallel to the plane of sources (here, the xy-plane). The object moves through the system on a conveyor belt in the z-direction. This results in an effective trajectory of the source respective to the object, which can be described as a function z : [−π, (2t − 1)π ] → I z ⊂ R. Here t is the winding number of the trajectory i.e. the number of activations of one physical source λ within one period of z. In the following the cone beam rays are parametrized by the Cartesian coordinates (u, v) of their intersection with a virtual detector plane D(λ), containing the z-axis, perpendicular to the xy-plane and facing the active source λ at the right angle.
Two sets of coordinates in a plane are used throughout the paper: the Cartesian coordinates (x, y) and the ray based coordinates (λ, u, l), where the source λ ∈ [−π, π ) along with the horizontal coordinate in the virtual detector plane, |u| u max , uniquely identify a ray in a plane, while the parameter l fixes a point along the ray (λ, u). The transformation from (λ, u, l) to the (x, y) coordinates and its inverse are respectively, where R denotes the radius of the ring of sources. The offset multi-source geometry over one effective trajectory period maps the spacial density function f (x, y, z) which vanishes outside of the bounded cylinder into the 3D ray transform g (λ, u, v) g(λ, u, v) = R 2 + u 2 + v 2 l 0 (u)+ l(u) l 0 (u)− l (u) dl f (X (λ, u, l), Y (λ, u, l), z(λ) + lv) (4) whereλ with v 1 (λ, u), v 2 (λ, u) projections of the first and last cylindrical detector row on the virtual planar detector D(λ) and u max = RR FOV / R 2 − R 2 FOV . The integration limits in (4) correspond to the intersection points of the planar ray (λ, u) with the transaxial field-of-view, where In the first paper in this series we defined the multi-sheet surface ζ 0 as an ordered set of surfaces (sheets) that are graphs of functions ζ s 0 , s = 1, . . . , S, where S = 2t is the number of sheets. The 2D fan beam transform on the multi-sheet surface ζ 0 (including all its sheets) was defined as with a 2D fan beam transform on a single sheet of the surface, ζ s 0 . Here f ζ s 0 (x, y) = f (x, y, ζ s 0 (x, y)), a 2D image selected from a function f defined on a 3D volume , as its values at the intersection of the volume and the sth sheet of the surface z = ζ s 0 (x, y). The rebinning of the offset multi-source geometry data to the multi-sheet surface ζ 0 is prescribed by the rebinning function V 0 : The rebinned data g 0 approximates the following mixed 2D fan beam transforms on the multi-sheet surface ζ 0 where i denotes the ith segment of the ray (λ, u). The superscript r = 1, . . . , t refers to the rth activation of the same physical source, λ ∈ [−π, π ), within one period of the effective trajectory,λ ∈ [−π, (2t − 1)π ), and is implicitly defined through the relatioñ The sheet ζ s 0 on which f is integrated along the ray segment i is determined by the mapping μ i (r) = s, which in the interval i , assigns the ith ray segment (λ, u, v, l), l ∈ i of the r(λ)th ray (λ, u, v) to the sheet ζ s 0 . The following implicit relation holds between the fan beam transforms on the S = 2t sheets and the t mixed fan beam transforms on the multi-sheet surface where (λ c , u c ) := (λ + π − 2 atan2(u, R), −u) 3 denotes the complementary ray to (λ, u) in a plane (reversed ray (λ, u)). As eachg (r) 0 (λ, u) andg (r) 0 (λ c , u c ) can be approximated using the rebinned offset multi-source data, (12) provides an approximation to p 0 .

Reconstruction on a multi-sheet surface
By definition the fan beam transform on the multi-sheet surface p 0 , λ ∈ [−π, π ), |u| u max is a sum of fan beam transforms on individual sheets p s 0 , s = 1, . . . , S. Substituting the expression for the transform on an individual sheet (7) into the definition of p 0 (6) and interchanging the order of summation and integration using linearity of those operators reveals that p 0 is also the fan beam transform of the superposition of the images on all sheets of the surface  Thus the task is to reconstruct the image f ζ 0 := S s=1 f ζ s 0 containing superposed information about the images on all individual sheets for subsequent deconvolution (unmixing the sheets). As in standard rebinning methods this stage is inherently parallelizable as the 2D problems on multi-sheet surfaces are independent from one another.

Sampling of the multi-sheet surface
In the ideal continuous case, the domain of the full scan fan beam transform in the plane for each ray contains the corresponding complementary ray and vice versa. Hence in the continuous case, it does not matter if the integral is taken along the primary or the complementary ray. However, as soon as we move to discretization, the exact complementary ray is in general not present in the domain of the discrete transform. A notable exception is what Natterer terms the PET geometry [15], in which the source and the detector positions coincide. The offset multi-source geometry could be in principle realized so that its projection on the transaxial plane is a PET geometry, by choosing the radius of the cylinder of detectors to be equal to the radius of the ring of sources and placing the sources and detectors at the same locations in the xy-plane but axially displaced from one another. The PET geometry is theoretically very appealing for its sampling efficiency [15], also it can be rearranged to parallel projection data. If the offset multi-source geometry is not a PET geometry, the assignment of ray segments to sheets (11), μ i , i = 1, . . . , S, in the mixed fan beam transform (9) will result in different effective sampling on different sheets. Figure 1 illustrates the sampling on sheets of a two-sheet surface.

Filtered backprojection on multi-sheet surface
As p 0 is a full scan fan beam transform, we can invert it using any full scan reconstruction formula.
As explained in the previous section after the discretization in generalg (r) 0 (λ c , u u ) is not available for a given (λ, u) thus it will need to be interpolated from neighboring values of g (r) 0 . The second problem is the cylindrical detector in the RTT geometry, in particular the non-uniform sampling it effects on p 0 (λ, u) along u. Mapping to arc instead of flat detector fan beam geometry also results in non-uniform sampling.
In contrast, in parallel geometry or PET geometry (which is a special case of parallel geometry) with even number of projections over 2π , for each ray (λ , u ) (λ ∈ [−π, π ) and |u | u max are the parallel projection angle and detector coordinate, respectively) the complementary ray (λ c , u c ) := (π +λ , −u ) is also contained in the domain of the discretized transform, thus the discretized sinogram preserves the symmetry of the Radon transform R f 2D (λ , u ) = R f 2D (π + λ , −u ). Therefore a preferable solution is to interpolate the parallel projection data or the PET data p 0 (λ , u ) and reconstruct using parallel filtered backprojection.
From (12) we obtain for the parallel transform on multi-sheet surface As the mixed fan beam transform is non-redundant, the interpolation ofg (r) 0 (λ c , u c ) and the fan to parallel interpolationg ,(r) 0 should use only proximate rays going in the same direction as the interpolated ray.
Due to the symmetry it is sufficient to reconstruct using only the first half of the parallel sinogram, λ ∈ [−π, 0), where p ,F 0 denotes the filtered parallel projection with h ramp the ramp filter and U (λ , x, y) = x cos λ + y sin λ .

Errors in the rebinned sinograms
In the previous section we discussed the reconstruction of the exact up to sampling full scan fan beam transform on the multi-sheet surface (12), p 0 . In practice however, we only have the rebinning approximation to p 0 . As the rebinned cone beam rays, g 0 , in general do not exactly lie in the surface (in contrast tog 0 ), this approximation suffers from errors in the axial direction. One consequence of the axial error is axial blurring in the reconstructed volumetric image. Incorporating blur model into the axial deconvolution problem significantly reduces this effect as explained in the next section. Another manifestation of the rebinning errors is inconsistency of the rebinned sinogram meaning that it is not a Radon transform of any nonnegative compactly supported function. Filtered backprojection of such an erroneous sinogram introduces an additional numerical error to the reconstruction on the multi-sheet surface. In contrast to the axial blur, these errors are non-local in the xy-plane.
The axial resolution phantom defined in section 8.3 was designed to reveal the axial errors due to rebinning. The error is expected to be the largest for a multi-sheet surface cutting through the phantom where the frequency of variation of the attenuation coefficient is the highest, here ζ 378 (378 refers to the index of multi-sheet surface in figures 4 and 5). Figure 2(a) shows the data rebinned to the multi-sheet surface ζ 378 , g RTT 378 , presented in a format we call an RTT sinogram with the abscissa indexing the full ring of sources, λ ∈ [−π, π ), and the ordinate the full ring of detectors, φ ∈ [−π, π ). The 2π -parallel sinogram, g 378 , obtained from g RTT 378 with a fan to parallel rebinning algorithm is depicted in figure 2(b). The rebinning approximation to p 378 in figure 2(c) was computed as the sum of rebinning approximations to the primary and complementary rays, g 378 (λ , u ) + g 378 (π + λ , −u ), q.v. (14). Due to the Radon transform symmetry, it is sufficient to approximate p 0 in the interval [−π, 0). The rebinning approximation to p 378 is an inconsistent sinogram, which gives rise to streak artifacts concentrated around the edges of the cylinders in the reconstructed image in figure 2(d).

Figures 3(c) and (d)
show the difference between the model of the blurred phantom image, (Af ) 378 (figure 3(a)) and the image reconstructed from the data rebinned to the multi-sheet surface ζ 378 (figure 2(d)) and their corresponding π -parallel sinograms, respectively. Thus figure 3(d) depicts the component of the error in the data not captured by the PSF model. Nonetheless, in the image space the densities seem to be quite faithfully reproduced with the PSF convolution model and the errors take the form of streak artifacts concentrated around the edges of the cylinders. Some image post-processing of the reconstructions on multisheet surfaces could be employed prior to the axial deconvolution to mollify those artifacts. Alternatively, algebraic reconstruction methods enhanced with e.g. sparsity constraints could be used for the reconstruction on the multi-sheet surfaces instead of filtered backprojection. Practical applications would however require extremely efficient implementation of the latter. In section 7.3 using local intra-block regularization during the axial deconvolution was suggested as a way of suppressing local variations in the xy-plane at a cost of solving  larger least-squares problems with a single right-hand side. Unfortunately, the size of the blocks directly relates to both the effectiveness of streak removal and the increase of the computational cost.

Idealized model of the convolution matrix
In section 3 using linearity of the fan beam transform we demonstrated that p ζ m is the fan beam transform of the superposition of images on all sheets of the multi-sheet surface ζ m . Therefore, we can first filter backproject the fan beam transform on the multi-sheet surface p ζ m (reconstructing an image f ζ m = S s=1 ζ s m which is a superposition of images on all the sheets of the surface) and subsequently in the image domain solve the system of simultaneous equations (17) to recover the volumetric density function f where M Z := {m : ζ m ∩ = ∅} is the set of indices of multi-sheet surfaces intersecting the bounded cylindrical volume . The equation system (17) is coupled due to intersection between surfaces filling the volume . A unique solution exists provided the multi-sheet surfaces sample the volume densely enough. Furthermore, the system (17) decouples into independent blocks for each fixed value of (x, y) as it is the case for the interpolation in single-sheet surface rebinning.
In the numerical procedure the function f is reconstructed on a voxel grid. We define a grid of cubical voxels, wherex = (x 1 , x 2 , . . . , x I ),ŷ = (y 1 , y 2 , . . . , y J ) andẑ = (z 1 , z 2 , . . . , z K ) are the centers of the voxels along x, y, and z axis, respectively. On this grid we have the discretized versions of the volumetric density functionf performs the rounding to the nearest integer and h x , h y , h z denote the voxel size along the respective axes.
Taking advantage of the independence of the equations between different pairs (x i , y j ), i ∈ 1, . . . , I, j ∈ 1, . . . , J, the discrete version of (17) decouples into I × J much smaller mutually independent systems of equations where Here χ b : Z → {0, 1} is the following set indicator function . . , S is the voxel valued multi-sheet surface. All of the independent systems (18) can be solved in parallel.
In practice we do not have the exact image on the multi-sheet surfaces fζ m , but only its approximation obtained from reconstructing the data rebinned to ζ m which results in erroneous right-hand sideb i j = b i j + i j and hence inconsistency of the system (18). In the next section we will show that better reconstructions can be obtained considering backprojection based models of the axial convolution matrix, which will result in a smaller data misfit.

Backprojection based model of the convolution matrix
The axial point spread function PSF m := PSF ζ m associated with the surface ζ m captures the distribution of the axial distances of the rebinned rays from the surface at every point (x, y). As the PSF m (x, y) depends only on z(λ), ζ m (x, y), V m (λ, u), whenever the problem of computing (ζ m (x, y), V m (λ, u)) can be reduced to computing a single function pair (ζ 0 (x, y), V 0 (λ, u)), the corresponding PSF m can be obtained from the PSF 0 by means of rotation Recall, that the reconstructed image on the multi-sheet surface is where B is the filter backprojection of the general form If we disregard the filtering (which corresponds to filtering with the Dirac delta kernel), we are left with backprojection, which is a local operation in the xy-plane in the sense that only those rays crossing the point (x, y) contribute to the backprojected image at this point each with the backprojection weight W (λ, x, y). We define the PSF 0 (x, y) as the ordered set of histograms PSF ζ s 0 (x, y) of the axial distances of the rebinned rays from the sheet ζ s 0 at the point (x, y) weighted by the corresponding backprojection weights where integrates only over the sheet ζ s 0 , i.e.λ : μ i λ +π 2π Notice, the use ofζ s 0 andV 0 which refer to the discretized version of ζ s 0 and quantized version of V 0 , respectively. By quantization we mean that V 0 has been projected on the cylindrical RTT detector, binned to the detector rings and subsequently projected back onto the virtual central detector plane D(λ). As PSF is used to setup a convolution matrix for the discrete problem, using the discretized quantities results in a model which closer matches the data.
This definition allows us to retain the benefit of independence of the deconvolution problems in the xy-plane. For each point (x i , y i ) in the grid we individually assemble the convolution matrices according to the following formula where Z ext ⊃ Z is an extension of the Z interval to accommodate the side lobes of PSF ζ m . Note, that (19) and (25) have the same form, thus formally the idealized convolution matrix can be thought of as constructed using (25) and (26) This allows for unified treatment of the idealized and PSF model convolution matrices.
Obviously, the PSF should be computed with backprojection formula matching the one used in the filtered backprojection on the multi-sheet surface. In the following we are going to give explicit expressions for the PSF models for different backprojection formulae.

Full scan fan beam backprojection
In section 3.2 of [3] we introduced a concept of a non-redundant fan beam transform. In a nutshell, the full scan fan beam transform in a plane is a sum of two non-redundant fan beam transforms, in which each ray is measured exactly once. The angular ranges of those transforms are complementary, [λ 1 , λ 2 ), [λ 2 , λ 1 + 2π ), respectively and they depend on the considered point (x, y). The backprojection formula for the non-redundant fan beam transform at (x, y) reads 1 2 U (λ, x, y)).
The PSF model for the fan beam backprojection is obtained from (24) substituting the nonredundant backprojection weights above PSF fan The trick of reparametrizing the continuous form saved us the cumbersome interpolation of the complementary rays in the discrete fan beam backprojection (see discussion in section 3.2). While discretization of (27) does not exactly match the discrete fan beam backprojection it turns out the PSFs are rather insensitive to the exact numerical formula and in any case due to the complementary rays interpolation, any PSF model derived for that backprojection would be approximate.

Full scan parallel beam backprojection
While constructing the PSF for the parallel beam backprojection we need to revert to approximation of the distance of the parallel ray to the surface at a given point, as in general the projections of the approximating cone beam rays on the xy-plane will not cross this point.
To approximate PSF ζ s 0 at a chosen point (x, y) we proceed as follows.
(i) For each of the anglesλ we find the corresponding detector u -coordinate of the parallel ray through (x, y), r (λ , u = U (λ , x, y)).
(iv) Approximate the distance of the parallel cone beam ray r cone, (λ , u ) to the sheet ζ s 0 at the point (x, y) as Using the above approximation we obtain the following expression for the PSF where integrates overλ : μ i λ +π 2π + 1 = s with i indicating the i segment of the parallel ray (λ , u ), which corresponds to the sheet ζ s 0 . In the above formula we chose to apply the indicator function to the distance of each of the approximating cone beam rays separately, as this again results in a model more closely resembling the data. Figure 5 shows the fit of the parallel backprojection PSF model and the reconstruction on the multi-sheet surfaces for the axial resolution phantom. Before constructing the convolution matrix the PSF ζ s 0 (x, y, :), s = 1, . . . , S have been smoothed by convolution with a rectangular window of length 3. The data fit has been greatly improved in comparison to the same experiment without the PSF model in figure 4. The PSF fan ζ 0 obtained with (27) is almost identical to the PSF ζ 0 , hence we did not include the corresponding plots of the data fit and henceforth we will restrict ourselves to the parallel backprojection.

Uniqueness of the deconvolution
In general (18) and so (31), are not square but rectangular, overdetermined systems of possibly inconsistent equations, and hence they have to be solved in the least-squares sense. In order to obtain a uniquely solvable least-squares problem, the system matrix has to be of a full column rank which in general holds for a sufficiently dense axial sampling with the multisheet surfaces. For the idealized convolution model (19), this can be intuitively quantified as that each of the voxels has to be intersected by at least S different sheets, which we believe in general is also sufficient for the PSF model case. For a sufficiently dense sampling then the condition numbers of the convolution matrices A i j will be bounded (and reasonably small) and the least-squares solution will approach the solution of the original problem as the norm of the error in the right-hand side i j 2 goes to 0. In practice in presence of the right-hand side errors it may be beneficial to quite generously oversample the problem to alleviate the errors.

Efficient implementation of the deconvolution
Apart from the obvious parallelism mentioned in section 5, there is a number of ways in which the solution of the systems (18) can be made numerically efficient.

Local data dependence
First we observe that the convolution matrices A i j are banded. Obviously, the band is larger for the convolution matrices assembled using PSF models. The band structure localizes the dependence of the solution on the right-hand side. Thus it is possible to split the large system obtained for the whole object into a series of small systems with matrices and right-hand sides corresponding to overlapping blocks of A i j , b i j , respectively. As a consequence, the deconvolution step can be implemented using a small buffer window at a given time holding only a fraction of the reconstructions on multi-sheet surfaces.

Assembling a computational block of the convolution matrix
For periodic effective trajectory, sampling of the volume with multi-sheet surfaces is also periodic. Due to the periodicity the banded convolution matrix decomposes into identical blocks (tiles), each corresponding to one period of the effective trajectory. Thus it suffices to assemble only one tile, then the large matrix A i j for the entire object can subsequently be built using shifted copies of this tile. We do not however want to build the matrix for the entire object but rather a smallest possible computational block for block-wise solution of the large system. In practice we may want to consider tiles corresponding to more than one effective trajectory period. This allows additional flexibility when choosing the block overlap.
First, we assemble the tile, A i j tile corresponding to the axial field-of-view of the scanner, Z tile , during t tile periods of the effective trajectory. To this end we use (25) but we consider only those surfaces,ζ m : m ∈ M Z tile that intersect Z tile : where Z ext.tile ⊃ Z tile extends Z tile to accommodate the side lobs of PSF ζ m . Figure 6 shows a fragment of A i j , A 2tile , corresponding to two consecutive A i j tile , stacked on the top of one another, while the second tile is also shifted by K tile = #Z tile .
We are going to use A 2tile to construct the matrix A for a block-wise solution of the large system. We choose a matrix [A i j block,0 | A i j block ] in such a way that it contains all the rows of A 2tile with non-zero entries in Z block (the part of the axial region Z, for which we want to solve with A) and crop all the zero columns, thus for [f i j block,0 |f i j block |f i j block,∞ ] T , the corresponding part of the solutionf i j , we have where b i j block is the part of the right-hand side of (18) with m ∈ M Z block . Solving the large system block by block, corresponds to solving for portions off i j which move through Z in K block size steps. In each stepf i j block,0 corresponds to the solution computed in the previous step, hencê f i j block,0 is already known and [f i j block |f i j block,∞ ] T can be obtained solving the rearranged system of equations f i j block,∞ corresponds to the part off i j , which is not completely determined by (30) (or (31)) and will be recomputed in the next step. In principle A i j block could be chosen differently: smaller or larger changingf i j block,∞ . Our choice corresponds to the smallest possible matrix capturing all the information available onf i j block . While solving (30) corresponds to almost independent solution of each of the blocks, the solutions of (31) are coupled between the blocks by the adjustment to the right-hand side. Therefore choosing to solve the system (31) instead of (30) carries the risk of error propagation between the blocks iff i j block,0 from the previous steps should be of low accuracy (e.g. a more challenging part of the image). The error inf i j block,0 will introduce an error to the right-hand side, block by block degrading the overall solution. In such a case it could be beneficial to solve the larger system (30), where the value off i j block,0 from previous step is not used. We would like to stress that when solving (30) we do not attempt to correct the solution from the previous step, just decouple it from the solution of the currently considered block. In the remainder of the paper we consider the system (30) for block-wise solution of the deconvolution problem, but all the results carry over to the system (31) without restrictions.

Using local similarity of convolution matrices
As the convolution matrices depend on slowly varying functions ζ m , V m and z, it is a valid simplification to evaluate ζ m and PSF ζ m on a coarser grid and compute the corresponding coarse grid matrices A. We then have to solve a much smaller number of systems of equations, each with many right-hand sides saving both memory and computational time.

Numerical solution of the deconvolution equation
Throughout this section we used a 32-threaded helix trajectory. The choice of trajectory influences the regularity of the sampling of the volume. An ideal sampling is achieved with rotationally symmetric surfaces which have rotationally symmetric PSF, [2]. In general, the more regular the sampling the better stability of the deconvolution step can be expected. This in turn allows us to use fewer multi-sheet surfaces per period of the effective trajectory resulting in smaller size convolution matrices saving storage and computational time. Figure 7 shows both the idealized and the PSF model convolution matrices evaluated on the coarse grid (5 × 5 pixels), for the blocks corresponding to the centers of the cylinders

Regularized L 2 solution
We consider a regularized version of the system (30) for a block-wise solution of the deconvolution problem. As mentioned in section 6.3 the convolution matrices used for blockwise solution in (30) can be evaluated on a coarser grid (of center points of k x × k y pixel blocks in the xy-plane). In this way the computational time can be reduced by solving k x k y times fewer systems each with k x k y right-hand sides Here, X and B are block matrices, which columns are solution and right-hand side vectors, respectively, corresponding to the pixels in the considered k x × k y block in the xy-plane, and · F denotes the Frobenius norm.
To take advantage of multiple right-hand sides it is necessary to factorize the system matrix. While QR factorization is more economical in storage, in the context of regularization a particularly useful factorization is the (G)SVD [8]. For each sparse convolution matrix, using a precomputed (G)SVD requires storage of two square dense matrices as well as the (generalized) singular values. Once the (G)SVD factorization is available we can compute for instance a Tikhonov regularized solution or a truncated (G)SVD solution for any value of the regularization parameter at a cost of order of dense matrix vector multiplication.
It is well known that imposing nonnegativity constraints when appropriate can result in a better quality reconstruction. In section 8.3 we used the nonnegative least squares algorithm by Lawson and Hanson [13] to solve the regularized version of (30). Unfortunately the method cannot benefit from multiple right-hand sides and moreover, even for a single right-hand side it converges rather slowly, which prohibits its use in practice.

Regularized L p /L q solution
Deconvolution problems are frequently tackled minimizing L p , 1 p < 2 norms In particular for q = 1 and L the forward difference operator we obtain the total variation (TV) regularization.
Minimization problem (33) can be solved e.g. using iteratively reweighted least-squares method (IRLS). The convergence of the scheme has been proven in [18] for 1 < p < 3, and p = 1 under uniqueness assumption on the solution. IRLS algorithm minimizes the residual norm r p p solving a series of least-squares problems, where in the kth step of IRLS the residual r k+1 has been reweighted r k+1 |r k | (p−2)/2 (multiplication is componentwise) using the residual from the previous step, r k . As reweighting changes the system matrix for each right-hand side individually, IRLS does not benefit from multiple right-hand sides in a general case. If however, we assume that the neighboring right-hand sides and so the solutions are sufficiently similar, we can use the same weights for the entire block.
The resulting block iteratively reweighted least-squares (BIRLS) method is summarized in algorithm 1. The BIRLS weights in the kth step, D k , are chosen as either where r i j k , η i j k are the columns of R k , η k , respectively, and w i j are normalized coefficients weighting each equation in the k x × k y block. Note, that (34) corresponds to using the constructed averaged residual to compute the BIRLS weights, while (35) to computing the IRLS weights for each individual right-hand side and then constructing the BIRLS weights as their weighted average.
The cost of BIRLS iteration is dominated by the cost of solution of the least-squares problem including the QR factorization (as opposed to the L 2 case where factorization could be precomputed once for all). In addition updating techniques could be used to further reduce the cost [17]. Our experiments have shown that only very few, around 3, BIRLS iterations are necessary to obtain a good solution.

3D intra-block regularization
Until now each of the right-hand sides in the k x × k y block in the xy-plane has been treated as independent and regularization was only applied in the axial direction. However, it may be beneficial to impose some regularity conditions on the solution locally in the xy-plane. Such conditions will couple the until now independentf i j block for different values of i, j. This approach is therefore not feasible for the entire transaxial plane as the resulting least-squares problems would be too large. However, within one k x × k y block we obtain a least squares problem with the matrix of k x k y times the size of the size of A in (32) where the regularization matrix comprises L z , L y , L x penalty operators applied along the axial and transaxial directions and x ∈ R k x k y n and b ∈ R k x k y m are vectors obtained by concatenation of the columns of X and B in (32), respectively. The 3D intra-block regularization can be seen as an alternative to image processing, e.g. filtering of images reconstructed on the multi-sheet surfaces prior to the deconvolution.

Reconstruction
In this section we show examples of volumetric images reconstructed with multi-sheet rebinning method. In our phantom experiments we used an offset multi-source geometry with a static ring containing 1152 sources and 10 rings each containing 1400 detectors. The axial offset of the first detector ring from the plane containing the sources was 10 mm and the axial spacing between the consecutive detector rings 2 mm. The radii of the source and detectors rings were 600 mm and 450 mm, respectively and the radius of the field-of-view R FOV was 400 mm. We assumed a constant axial translation of the scanner of 16 mm per period of the firing sequence, 2π . The effective trajectory was a 32-threaded helix. The images were reconstructed on an 1 × 1 × 0.5 mm grid. The simulated data was obtained with the Siddon algorithm using voxels of half the size of those used for the reconstruction. The detector aperture was emulated shooting rays to nine equally distributed points on the detector pixel. Subsequently 3% of Gaussian pseudorandom noise was added.
For the reconstruction of the phantom images, the simulated cone beam data was rebinned to 96 multi-sheet surfaces per effective trajectory period (corresponds to three times the number of the transaxial slices through the volume we want to recover). The reconstruction on the multi-sheet surfaces was performed using parallel filtered backprojection following unidirectional fan to parallel rebinning with 800 parallel rays and 1258 angles over 360 • range. The PSF model (25), (28) used for construction of the convolution matrices has been evaluated on a grid of 5 × 5 pixel blocks in the transaxial plane. The PSFs have been posteriori smoothed with a rectangular window of width 3. For all forms of the deconvolution problem we used the regularization parameter calculated in section 8.1, but it was multiplied with an image specific factor for the TV regularization. The L 1 and TV regularized solutions have been obtained with at most three BIRLS iterations subject to relative residual change of more than 10 −3 .

Regularization parameter
Underlying our heuristic choice of the regularization parameter for the deconvolution problem (32) is an assumption of spatial variation of the standard deviation of the data in the xy-plane, proportional to the axial deviation of the rays from the multi-sheet surfaces. The axial deviation from the multi-sheet surface is given by function where Q q has been discussed in detail in [3]. As the convolution matrix corresponds to multi-sheet surfaces at different rotations (throughout the range [−π, π )), Q q is subsequently radially averaged over annuli in k 2 x + k 2 y radial steps, and normalized to have average value of 1 (scaling such thatQ q on the annuli containing the radius R FOV / √ 2 becomes 1). Finally, we choose the spatially varying standard deviation as σQ 2 (x, y) (note that its average value is σ , the estimated standard deviation of the error between the PSF model and the reconstruction of the phantom on multi-sheet surfaces). Under assumption of zero mean and a diagonal covariance matrix, (σQ 2 (x i , y j )) 2 I, the expected value of the error norm Ax−b 2 is σQ 2 (x i , y j ) √ m(i, j), where m(i, j) is the number of rows of the convolution matrix A, (30). The spatially dependent value of the regularization parameter has been evaluated for the axial resolution phantom defined in section 8.3 and is shown in figure 9. It was then used for all reconstructions from simulated and real data for all choices of norms and penalties.
In particular in the L 2 case it is possible to evaluate the regularization parameter on fly using criteria such as L-curve to adapt to the local content of the image. In the multiple right-hand side case this could be done either for each column of X separately resulting is a different regularization parameter for each column, or for all the columns simultaneously resulting in one parameter for the entire block. While our experiments indicate that automated evaluation of the optimal regularization parameter can be unreliable, it can be part of the post-processing step, where the user specifies the part of the image to be deconvolved with different values of the regularization parameter. Note that due to the locality of the deconvolution problem in all three dimensions of space, any desired part of the image can be recomputed for any regularization parameter value without necessity of recomputing the entire volume.

Impact of the regularization on the discontinuities of the solution
The effect of regularization on the jump singularities can be evaluated as follows. We consider Heaviside function, H(z), as a model jump singularity. Using the deconvolution matrix A, we construct the ideal data for this model, b = AH. Subsequently, we solve the regularized deconvolution problem, and compare the solution to the Heaviside function. The same approach has been invoked in [14] to investigate the impact of window parameter in the windowed singular value solution of the deconvolution problem for functions with point singularities. Figures 10(a) and (b) show the L 2 and TV regularized solutions for deconvolution problems (30) with A 400,400 , A 512,288 and A 632,168 and L the forward difference operator. The regularization parameter α for L 2 is shown in figure 9(b) and 0.1α is used for the TV regularization. TV regularization recovers sharper jumps without over or under shooting in contrast to the L 2 regularization. The reconstructed jumps are the sharpest in the center of the transaxial plane, and diffuse away from the center.

Axial resolution phantom
As all rebinning methods, multi-sheet surface rebinning suffers diminished axial resolution due to the rebinning errors. Our first phantom is therefore tailored to demonstrate the limitations of the axial resolution of the method. The axial resolution phantom contains three cylinders, which centers have been equidistributed along the radius of the transaxial field-of-view of the scanner, see figure 11(a). The attenuation coefficient of each of those cylinders has been chosen to vary in decreasing frequency along the axial direction from every 2 to 5 mm finishing with a longer cylinder of constant attenuation, see figure 11(b). Placing the oscillating cylinders at different positions, tests the axial resolution throughout the transaxial plane. The best resolution is expected in the central region as the bands of the corresponding convolution matrices are the slimmest, degrading outwards due to the increasing bandwidth. This phantom has been used to illustrate the errors in the rebinned sinograms and their reconstructions in section 4 as well as to evaluate the data fit of the convolution matrix models in section 5. The reconstructions of the volume with different strategies of solving the deconvolution problem are illustrated in figures 12-14. Contrast in all reconstructed images have been set equal for better comparability. Figure 12 shows the superposition of the reconstruction on xz-planes crossing through the centers of each of the cylinders, y = 400, 288, 168. The profiles of the solutions corresponding to cross-section through the centers of each of the cylinders are shown in figure 14 for a selection of solution methods. A selected transaxial slice through the volumes in the high oscillatory region at z = 136 is depicted in figure 13.
Comparison of figures 12(a) and (b) immediately confirms that using the PSF model yields superior results. With the idealized deconvolution it is not possible to resolve the 3 mm spaced discs of the middle cylinder and not even the 5 mm spaced ones of the outward cylinder, while the 3 mm spaced discs in both cylinders can be made out for the PSF model. Minimizing the L 1 norm provides slightly sharper image but is more prone to overshooting in some higher frequency parts of solution as seen for the outward cylinder in figure 14(c). While the streak artifacts are preserved by the L 1 reconstruction (figure 13(c)), they are largely suppressed by the TV solution (figure 13(e), regularization parameter 0.1α) which provides sharpness similar to L 1 , but seems slightly more stable, see figure 14(c). The 3D intra-block regularization in figure 13(d) slightly smoothed the streak artifacts but the intensity variation of the background remained, while the nonnegativity constraint reduced the streaks considerably, see figure 13(f).

Clock phantom
Our second example is the sphere clock phantom [20]  TV regularization effectively eliminates the oscillations in the L 2 solution as well as the streak artifacts.

Real data
Finally, we would like to demonstrate the performance of the MSSR on the real data. The data were obtained with the Rapiscan Systems RTT80 machine. Figure 17 shows the computed x-ray images (integrated along the thinnest dimension) of an L 2 and TV regularized reconstruction of a five string bass guitar.

Conclusions
We presented a new family of methods for reconstruction from axially asymmetrically truncated cone beam projections. The methods combine analytical and numerical ideas, exploiting the linearity of ray transforms to reconstruct data on multi-sheet surfaces and subsequently perform axial deconvolution to recover the volumetric image from a stack of images reconstructed on multi-sheet surfaces. In the first paper in this series we discussed the approximation of the data on multi-sheet surfaces. We set up a variational problem for obtaining an optimal multi-sheet surface and rebinning function and showed that it can be solved by a globally convergent iteration in the case of quadratic cost function. In this second paper in the series we have shown how the volumetric image can be reconstructed from data rebinned to multi-sheet surfaces. We discussed in detail the reconstruction on a multisheet surface and the effect of the rebinning errors. We set up a deconvolution problem for the recovery of the volumetric image from a stack of the images on multi-sheet surfaces.
We developed backprojection based models of axial convolution matrix, which substantially reduced the data model misfit and hence improved the quality of the solution. We outlined an efficient implementation of the deconvolution step and explored different strategies for its solution. Finally, we demonstrated the performance of the method on both the simulated data and real data collected with a Rapiscan RTT80 cone beam scanner. An interesting future research direction is combining the ideas of rebinning with sparsity enhanced reconstruction from subsampled projections. This would allow for even faster data acquisition provided the imaged objects have sparse representation in the chosen basis. Furthermore, fitting fewer projections to a surface potentially reduces the rebinning errors.