Multi-sheet surface rebinning methods for reconstruction from asymmetrically truncated cone beam projections: I. Approximation and optimality

The mechanical motion of the gantry in conventional cone beam CT scanners restricts the speed of data acquisition in applications with near real time requirements. A possible resolution of this problem is to replace the moving source detector assembly with static parts that are electronically activated. An example of such a system is the Rapiscan Systems RTT80 real time tomography scanner, with a static ring of sources and axially offset static cylinder of detectors. A consequence of such a design is asymmetrical axial truncation of the cone beam projections resulting, in the sense of integral geometry, in severely incomplete data. In particular we collect data only in a fraction of the Tam–Danielsson window, hence the standard cone beam reconstruction techniques do not apply. In this work we propose a family of multi-sheet surface rebinning methods for reconstruction from such truncated projections. The proposed methods combine analytical and numerical ideas utilizing linearity of the ray transform to reconstruct data on multi-sheet surfaces, from which the volumetric image is obtained through deconvolution. In this first paper in the series, we discuss the rebinning to multi-sheet surfaces. In particular we concentrate on the underlying transforms on multi-sheet surfaces and their approximation with data collected by offset multi-source scanning geometries like the RTT. The optimal multi-sheet surface and the corresponding rebinning function are found as a solution of a variational problem. In the case of the quadratic objective, the variational problem for the optimal rebinning pair can be solved by a globally convergent iteration. Examples of optimal rebinning pairs are computed for different trajectories. We formulate the axial deconvolution problem for the recovery of the volumetric image from the reconstructions on multi-sheet surfaces. Efficient and stable solution of the deconvolution problem is the subject of the second paper in this series (Betcke and Lionheart 2013 Inverse Problems 29 115004).


Introduction
There are obvious situations in which fast acquisition of the cone beam computerized tomography data is of utmost importance. In medical applications for instance, cardiac or lung imaging at present require gating techniques to keep up with the rapid motion of the organ. However, gating diminishes the quality of the reconstructed images potentially compromising their diagnostic value. Another time critical application is baggage security screening where the scanning speed directly relates to the equipment costs for the airports and clearance time for air passengers. Fast data acquisition also enables study and control of processes occurring at short time scales, such as the flow of oil-air-gas mixtures [23].
A conventional state of the art x-ray cone beam system uses one source of radiation positioned opposite an array of detectors. The detector and source assembly is then moved mechanically relative to the object to be imaged: typically the source and detector rotate around the object while the object is translated resulting in a helical trajectory relative to the object (HCB). The rate at which tomographic images can be acquired by such a system is limited by the rate of rotation of the assembly supporting the source and detector array. One way to speed up data acquisition is to replace the mechanically rotating gantry by a stationary ring of sources, which can be quickly switched electronically, and multiple stationary rings of detectors. Switching the sources on and off in a sequence gives the same effect as a rotating source, though they can also be activated in any other order. In such design however, no ring of detectors can be placed in the same plane as the ring of sources as it would obstruct the beam. This means that attenuation along rays that make less than a certain limiting angle to the plane of the sources cannot be measured. In fact, the axial offset of the detector from the plane of sources can be substantial (of the order of the detector axial length). We call such a geometry an offset multi-source geometry [15,16]. The offset multi-source geometry data are nontrivially axially truncated and hence they require new reconstruction algorithms different from those devised for the standard cone beam CT. Potentially, the detectors could be placed on both sides of the ring of sources, resulting in an axially symmetric geometry still missing the rays within the limiting angle. In this work we concentrate on the more challenging case with the detectors on only one side of the ring of sources.
The basis idea of rebinning methods [4-6, 8, 10-12, 14, 17, 24, 25] is to reduce the original three-dimensional (3D) reconstruction problem to a series of approximated twodimensional (2D) reconstruction problems e.g. corresponding to data on transaxial slices from which the volumetric image is subsequently interpolated. Both the 2D reconstruction and the interpolation can be very efficiently implemented (e.g. in hardware). Yet another benefit of rebinning methods is their local data dependence enabling on-fly reconstruction i.e. without necessity to acquire the measurements of an entire object before initiating the reconstruction. Those qualities make rebinning methods particularly attractive for applications with real time requirements such as security screening.
In this work we introduce a new family of rebinning methods, multi-sheet surface rebinning (MSSR) methods [1], for reconstruction of offset multi-source cone beam CT data. MSSR methods approximate and reconstruct data on multi-sheet surfaces and use the linearity of the ray transform to obtain the volumetric image through axial deconvolution. They retain the efficiency and local data dependence of the rebinning methods while allowing for high resolution reconstruction from, in the sense of integral geometry, severely limited data. A version of MSSR methods using parametric surface is presented in [2]. The choice of rotationally symmetric double-cone surface allows to derive a John's equation based correction for the rebinned data in the spirit of [9,13]. The correction results in effectively rotationally symmetric system. As a consequence, the corresponding deconvolution matrices have Toeplitz like structure and vary only with the distance to origin in the transaxial field-of-view, allowing very efficient and robust deconvolution.
The remainder of this paper is organized as follows. In the next section we develop a formal description of an offset multi-source geometry. In section 3 we introduce the concept of non-redundant fan beam transform, which will play an essential role in the derivation of MSSR methods. In section 4 we briefly review single surface rebinning methods for standard helical cone beam geometry and their generalization to offset multi-source geometry developed in [4]. The rebinning to multi-sheet surfaces is introduced in section 5. In section 6 we discuss an optimization problem for the choice of the multi-sheet rebinning surface and function pair and describe an alternating iteration for its solution, which convergence is proven in the appendix. We proceed to give some examples of two-sheet rebinning surface and function pairs in section 7. The axial deconvolution problem for the recovery of the volumetric image from reconstructed images on multi-sheet surfaces is derived in section 8. The detailed description of how the deconvolution can be efficiently and stably performed is subject of the second paper in this series [3]. We conclude with a summary of the results and a preview on the second paper. Numerical reconstructions are deliberately postponed to the second paper.

Offset multi-source geometry
Figure 1(a) shows the offset multi-source geometry as deployed in the Rapiscan RTT80 scanner [15,16]. As opposed to the conventional cone beam geometry the RTT geometry features a static ring of sources, and multiple static rings of detectors. The grey area on the detector cylinder indicates the part of the detector illuminated by the active source and the axial size and offset of the detector dictates the axial truncation of the x-ray cone. In the RTT geometry the ring of sources and the rings of detectors are each located in a plane parallel to the xy-plane, while the object moves through the system on a conveyor belt in the direction of z-axis.

Effective trajectory
Formal description of the offset multi-source geometry requires introduction of notion of effective trajectory to account for the freedom of choice of the firing order of sources in the ring.
Letλ j ∈ˆ := atan2 s j y , s j x : j ∈ 1, . . . , N s ⊂ [−π, π ) 3 denote the unique angular coordinate of the source in a ring, where triple s j x , s j y , s j z denotes the x, y and z coordinates of the source and the N s the number of sources in the ring. Without loss of generality we assume that the sources are located in the xy-plane, i.e. s j z = 0. We assume that each of the sources in the ring is activated exactly t times within tN s long source firing sequence,ˆ A = (λ i ),λ i ∈ˆ , and thereafterˆ A repeats itself, which results in a tN s periodic scanning process. In order to capture the translation of the system with respect to the object at each source activation, we introduce the axial displacement functionẑ :ˆ A → R.
We define the effective trajectory as a continuum generalization of the axial displacement functionẑ(λ i ),λ ∈ˆ A . To this end we identify the sources in the firing sequence with angular parameter over multiple turns, where # denotes the cardinality of the set. Asˆ A is periodic with the period tN s , for its continuum equivalent it is sufficient to consider λ ∈ I λ := [−π, (2t − 1)π ] and we can regard t as a winding number of the trajectory. We now define the axial displacement function z acting on an interval, z : I λ → I z ⊂ R, and we require that this function at the pointsλ i assumes the values z(λ i ) =ẑ(λ i ). The axial displacement function z is chosen from set of bounded variation functions that satisfy where the derivative is defined in the distributional sense. The constraint equalizing the total variation (TV) of the discrete and bounded variation functions eliminates unreasonable trajectories. Obviously, there are still infinitely many functions z satisfying these conditions. To eliminate the remaining ambiguity we assume z to be a linear interpolant ofẑ everywhere but on the intervals containing the jump points of z, where we use the assumption of constant scanning speed to define the interpolant. The so defined displacement function z fully captures motion of the scanner relative to the imaged object and it generalizes the concept of the trajectory to the system with multiple sources.

Ray parametrization
Although naturally the detectors are arranged on the cylinder for the sake of notational convenience we are going to parametrize measured rays by 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, figure 1(b). A pixel on the virtual detector is described by the Cartesian coordinates (u, v), where the orthonormal basis vectors 1 u and 1 v given by with the origin (u, v) = (0, 0) the orthogonal projection of the active source λ on D(λ). This parametrization is no restriction as the cylindrical detector can be easily mapped to the flat detector and vice versa.

Plane parametrization
We use two sets of coordinates in a plane for the points in the transaxial field-of-view 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), see figure 3. The transformation from (λ, u, l) to the Cartesian coordinates is described by the following formulae X (λ, u, l) = R cos λ + l(−R cos λ − u sin λ), Y (λ, u, l) = R sin λ + l(−R sin λ + u cos λ) (2) and the inverse transformation by where R is the radius of the ring of sources.

Asymmetrically axially truncated cone beam projections
In x-ray CT one seeks to reconstruct linear attenuation coefficient at every point in the object from a set of line measurements. The linear attenuation coefficient can be thought of as a function f (x, y, z) which vanishes outside of a bounded cylindrical field-of-view of the scanner Figure 2 shows the projections of the cylindrical detector on the virtual detector plane for HCB (a) and the RTT (b) geometries. In both cases, the area between the projections of The symmetry of HCB detector implies that the detector is practically unconstrained. Simply by adjusting scanner parameters such as slowing down the belt speed or equivalently symmetrically axially extending the detector, the entire Tam-Danielsson window [7,21] can be fitted in the active detector region. In contrast for the offset multi-source geometry, the lower bound v 1 is not altered in this way presenting a genuine constraint. While the upper bound v 2 could be eliminated in a way described above, we chose nonetheless to impose it to remain within the existing scanner parameters. The resulting detector constraints are for λ ∈ I λ , u ∈ [−u max , u max ]. The horizontal extend of the detector is assumed to be large enough to accommodate the entire field-of-view, with u max = RR FOV / R 2 − R 2 FOV . As a result of the detector constraints (5) the offset multi-source geometry over a period of effective trajectory z measures a set of line integrals with The integration limits in (6) correspond to the intersection points of the ray (λ, u) with the field-of-view 2D , where and l 0 (u) corresponds to the point of the closest approach of the ray (λ, u) and l(u) corresponds to the in plane distance along the ray (λ, u) from the point of the closest approach to the boundary of the field-of-view. Note that the axial detector coordinate v only appears in the scaling factor and in the last argument (the axial coordinate) of the spacial density function in (6), which is a consequence of using the cone beam data rather than the fan beam data. This separation of the axial and transaxial directions enables the rebinning methods discussed in section 5.

Complementary rays
In what follows we will need the concept of complementary rays. Any directional ray in a plane is uniquely described by its coordinates (λ, u), where λ ∈ [−π, π ) is the angular parameter of the vertex and u ∈ [−u max , u max ] is the coordinate on the virtual detector line. We call a ray (λ c , u c ) complementary to the ray (λ, u) if The complementary rays are simply the reversed rays, i.e. with the source and detectors roles swapped, see figure 4. The relation is symmetric thus (λ, u) is also complementary to (λ c , u c ).
Observing that 2l 0 is proportional to the length of the ray from the source to the complementary source, a simple geometrical argument yields which completes the complementary set of variables describing the point in a plane In the context of 3D cone beam geometry, we say that all rays shot within one period of the effective trajectory, for which orthogonal projection on the xy-plane is either (λ, u) or (λ c , u c ), are complementary. For the full 2tπ range scan, λ ∈ [−π, (2t − 1)π ), there are exactly 2t such rays.

Non-redundant fan beam transform
There are two conventionally used fan beam formulae corresponding to full scan (sources in the angular range [−π, π )) and short scan data (sources in the angular range [−π, 2δ), δ = arcsin(R FOV /R)). Both those formulae produce redundant data, as the integral along each line is taken twice in case of the full scan and at least once in case of the short scan fan beam transform. In [18] it was shown that for reconstruction of only a part of the field-of-view even fewer data are needed. Precisely, the part of the object which lies in the convex hull of the scan can be reconstructed from non-truncated fan beam projections.
In the spirit of [18], to reconstruct a single point P in the region of interest we need a set of non-truncated projections, each of which contains exactly one line passing through P, such that we have a (non-directional) line through P from every direction. As it can be seen in figure 4, there are infinitely many choices of a scan ranges which fulfil this condition e.g. [λ 1 , λ 2 ) or [λ, λ c ), however there is a unique minimal range. Namely, the one which corresponds to the end points of the line segment with the centre P, P being the point of the closest approach of the line to origin, [λ 1 , λ 2 ). Moreover, the maximal range corresponds to the complement of the minimal range, [λ 2 , λ 1 + 2π ). The length of the range for a point (r, θ ) depends only on its distance to the origin, r, while the range interval rotates with the polar angle θ , [λ 1 , λ 2 ) = θ − arccos r R , θ + arccos r R . We observe that the origin, r = 0, requires the largest minimal scan range of π to be reconstructed. On the other hand, the origin is the point with the smallest maximal range of π . We construct a non-redundant fan beam transform in which integral along each line is taken exactly once in the following way. For every point in the field-of-view we select only those rays from the full scan, which originate from the minimal angular scan range for this particular point. Though for each individual point sources only within the angular range of at most half scan (i.e. π ) are used, for the entire field-of-view the non-redundant fan beam transform requires sources from the angular range of 2π . Analogously, the transform obtained through selecting those rays originating from the maximal angular scan range for each point is a non-redundant fan beam transform, and moreover it is complementary to the minimal transform in the sense that the union of both transforms results in the full scan fan beam transform of angular range 2π . Obviously, none of the non-redundant transforms can be measured for every point in the region of interest by a fan beam scanner independently only their union. Assuming that we found a way to disentangle the non-redundant transforms, we ask what would be the appropriate backprojection weights for such a non-redundant fan beam transform.
Let f 2D (x, y), (x, y) ∈ 2D denote a 2D image. The full scan fan beam backprojection readsf where g 2D is the corresponding 2D fan beam data. We now decompose the full scan backprojection integral into the parts corresponding to the two non-redundant fan beam transforms where min := [λ 1 , λ 2 ) and max := [λ 2 , λ 1 +2π ) are the point P = (x, y) dependent minimal and maximal ranges, respectively. By the definition, for each ray from λ ∈ max through P, there is a complementary ray going in the opposite direction from λ c ∈ min through P. Thus we can reparametrize the second integral using the complementary set of variables (8) where we used shorthand notation λ c (λ, x, y) Here (14) is easily verified and the last equality (15) follows from (9).
x, y)), we conclude that the weights for the backprojection of the non-redundant fan beam transform with the minimal range arẽ Obviously, due to the reciprocity of the relation of being the complementary ray,f 2D (x, y) can be obtained backprojecting the maximal range non-redundant fan beam transform, using the same formula as in (16) but where the integration is performed over λ ∈ [λ 2 , λ 1 + 2π ]. The non-redundant backprojection formula lends itself to an easy interpretation. The weight for each ray integral is a sum of the fan beam backprojection weights corresponding to the ray itself and to its complementary ray. The additional weight compensates for the fact that each ray is backprojected only once.

Single-sheet surface rebinning methods
The basic idea behind rebinning methods is to rearrange (rebin) the set of cone beam data into possibly overlapping subsets, each approximating a complete set of data corresponding to some 2D problem. In this way the difficult 3D problem is reduced to a series of easier 2D problems, each one of which can be independently reconstructed using 2D filtered backprojection. As we start from the cone beam data, it is natural to approximate either the short scan or the full 2π scan fan beam data. In the last stage, the volumetric image is axially interpolated from the reconstructed 2D problems.
Rebinning methods seek an approximation which minimizes some measure of axial deviation of the rebinned cone beam rays from the desired fan beam rays, while maintaining the exact transaxial resolution. At the heart of this procedure lies the question how to choose the 2D problems, to be approximated. The natural way to approach this problem is to find a plausible surface, for which a complete data set can be found with a minimal approximation error. The first simplest idea was to rebin the cone beam data to transaxial slices, so-called single slice rebinning (SSR) [5,17]. An improvement was achieved using tilted planes tailored to the helical trajectory [11,12,14,24] and further an iteratively constructed nonplanar surface [6,10] and nonplanar parametric surface [25]. In [8] Defrise, Noo and Kudo presented a unified derivation of surface rebinning methods, posing the problem in variational framework, as minimization of square axial deviation functional with respect to the rebinning surface and function. They identified the previously introduced rebinning methods as special cases of this variational problem, using different parametrizations of the rebinning surface.
In [4] we generalized the approach in [8] to the offset multi-source geometry, by formulating a constrained variational problem, with constraints (5) on the rebinning function V . The interpolation was performed by specially tailored axial deviation informed scattered data interpolation algorithm. A major drawback of the constrained optimal surface rebinning is the nonzero curvature of the optimal surface ζ , which results in relative large axial deviation of the rebinned rays from the surface.

Rebinning equation for a (single-sheet) surface
For a given rebinning surface ζ λ 0 , let be a 2D image selected from a density function f defined on a 3D volume , as its values at the intersection of the volume and the surface z = ζ λ 0 (x, y). In the simplest case, when ζ λ 0 is a transaxial plane, the 2D image corresponds to transaxial slice through the volumetric function f .
Taking fan beam transform of the 2D image f λ 0 is then equivalent to the following fan beam transform on the surface ζ λ 0 where λ = π for the full 2π scan range and λ = π/2 + δ, for the short scan range, with δ = arcsin(R FOV /R) the fan beam half aperture. Following [8] we call λ 0 a rebinning centre of ζ λ 0 as for helical trajectory, λ 0 , is the centre of the helical segment, h[λ 0 − λ, λ 0 + λ], used for rebinning to the surface ζ λ 0 , where 2π h is the pitch of the helix. Since λ 0 is unique to the surface ζ λ 0 , it can serve as a surface index. Accordingly, the 2D image f ζ λ 0 can be reconstructed from p λ 0 using fan beam filtered backprojection. However, the surface data (17) cannot be measured. Instead, rebinning (18) is used to approximate the data p λ 0 on the surface ζ 0 with the measured cone beam rays. To avoid loss of transaxial resolution, rebinning methods use data along rays which have the same orthogonal projection on the xy-plane i.e. p λ 0 (λ, u) is approximated by The square root factor rescales the 3D cone beam data to consistent 2D fan beam data and ensures that rebinning for objects with density independent of z, f (x, y, z) = f (x, y, 0), is exact.
To each generalized ray on the surface ζ λ 0 , the rebinning function V λ 0 assigns an approximating cone beam ray, i.e. V λ 0 (λ, u) is the vertical coordinate on the virtual detector plane D(λ) of the cone beam ray (λ, u, V λ 0 (λ, u)) rebinned to the generalized ray (λ, u) on the surface ζ λ 0 .

Multi-sheet surface rebinning
The methods described in this paper improve the quality of the axial approximation on the constrained optimal surface rebinning methods [4] using a multi-sheet surface instead of a single-sheet surface.

Multi-sheet surface
We define a multi-sheet surface ζ as an ordered set of surfaces (sheets) that are graphs of functions ζ s , s ∈ 1, . . . , S, where S is the number of sheets. Obviously, a simple surface in the light of this definition is a single-sheet surface, a special case of a multi-sheet surface.
As in the single-sheet surface case we associate the multi-sheet rebinning surface ζ λ 0 with its rebinning centre λ 0 , which is a centre of the axial range of the trajectory segment used for rebinning to ζ λ 0 . For the helix it was the centre of the angular interval, but this is no longer the case for general trajectories.
For single-sheet surface and a short scan, choosing rebinning centres at equispaced angles guarantees uniform usage of cone beam projections from all angles. However, when a full scan angular range or its multiple (2tπ ) is used, uniform projection usage follows automatically and other selection criteria become important, such as uniform sampling of the volume by the associated rebinning surfaces. To this end, we are going to select a subsequence 0 ⊂ I λ of N c rebinning centres per firing sequence period, such that z(λ m ), λ m ∈ 0 are uniformly distributed along the z-axis. Evidently, 0 can be periodically extended with the period N c . Equation (19) then suggests uniform distribution of the surfaces in the volume. Clearly the less rotationally symmetric the surfaces are, the less uniform the effective sampling of the volume (e.g. constrained optimal single-sheet surfaces). We come back to this problem in section 7. Under the assumption of constant scanning speed such a choice of 0 amounts to selecting sources fired in equal time intervals.
In contrast to the single-sheet surface rebinning, where short scan can be used, the two-sheet surface rebinning requires cone beam projections from sources within the angular range of 2π , i.e. projections from all physical sources. In general the multi-sheet rebinning with S = 2t, t ∈ N sheets, will require data from angular range 2tπ . For rebinning to ζ λ 0 we consider projections acquired from sources in angular range (λ 0 −2tπ, λ 0 +2tπ ) active within the shortest possible axial distance from z(λ 0 ). For constant scanning speed this corresponds to using data from within one period of the firing sequence and to axial translation of the scanner [z(λ 0 ) − z t /2, z(λ 0 ) + z t /2), where z t = z((2t − 1)π ) − z(−π ) is the translation within one effective trajectory period. Henceforth, we refer to the set of sources rebinned to ζ λ 0 as As in the case of single-sheet surface, using other rays from outside of this shortest interval can result in improved signal to noise ratio [5,20], but we do not follow this line of thought here.

Symmetries of the rebinning surface and function
In general for the reconstruction of a function supported on a bounded cylindrical volume we need all the pairs (ζ m , V m ) := (ζ λ m , V λ m ), for which intersection with is non-empty, ζ m ∩ = ∅. Due to periodicity of the effective trajectory it is sufficient to consider the pairs (ζ m , V m ) within one period, m ∈ 0 . Then the remaining pairs, (ζ km , V km ) can then be obtained by the periodic extension As first observed by Heuscher [10] the invariance of helical trajectory under a combined rotation by θ and axial translation by hθ , where 2π h is the pitch of the helix, gives rise to symmetries which carry over to the rebinning pair, reducing the problem to finding a single rebinning surface and rebinning function (ζ 0 , V 0 ). This symmetry was used in [8] for simplifying the problem of finding an optimal surface and rebinning function for a standard helical cone beam geometry.
In fact the same invariance holds for an n-threaded helix, making it an appealing choice of firing sequence for an offset multi-source geometry. Heuscher's observation can be generalized to offset multi-source geometry arguing analogously to [8]. Let θ be a geometric transformation combining a rotation around z-axis by angle θ with an axial translation by z θ := z t θ/(2tπ ). The effects of applying θ in image and data space are, respectively This correspondence carries over to the optimal rebinning pair, whenever the scanning geometry satisfies the required invariance. Thus for every effective trajectory z(λ) invariant under θ , θ ∈ [−π, (2t − 1)π ), and for any detector constraints v 1 (λ, u), v 2 (λ, u) invariant under rotation by θ (such as the RTT geometry with n-threaded helix effective trajectory) for each sheet ζ s m of an optimal rebinning surface ζ m and the rebinning function reducing the problem to finding a single rebinning pair (ζ 0 , V 0 ) := (ζ λ 0 , V λ 0 ). We remark that the RTT geometry with two detector cylinders placed symmetrically with respect to the source plane also satisfies the required invariance. Exploiting the symmetry allows more memory efficient implementation, which is important if memory is scarce or accessing large data sets cannot be efficiently done as for example in graphic processing units. Henceforth, we are going to assume the θ symmetry but the extension to the general case is immediately obvious. In contrast to the standard HCB geometry, in the offset multi-source geometry due to the presence of positive constrains v 1 , v 2 > 0, we cannot expect the antisymmetry in y, ζ s 0 (x, y) = −ζ s 0 (x, −y), λ 0 = 0, to hold for a sheet of the rebinning surface (neither singlesheet nor multi-sheet) and consequently in general ζ 0 (R, 0) = 0. However, for a particular effective trajectory, the shape of the rebinning pair is preserved under axial and transaxial dilations.

Rebinning equation for the two-sheet surface
The simplest realization of the multi-sheet rebinning is rebinning to a surface ζ 0 with S = 2 sheets, ζ s 0 , s ∈ {1, 2}. In the following we frequently use labels: b (bottom) for s = 1, and t (top) for s = 2 to improve readability. In the two-sheet case we assume the winding number of the trajectory t = 1, i.e. every source inˆ is fired within one trajectory period exactly once before the pattern repeats itself. Consequently, λ ∈ [−π, π ), |u| u max throughout the section. On each sheet of the surface ζ s 0 we define the 2D fan beam transform, exactly as in the case of single-sheet surface (17) Furthermore, we define a 2D fan beam transform on the entire multi-sheet surface ζ 0 (including all its sheets) as a superposition of the fan beam transforms on all the individual sheets Figure 5(a) shows p b 0 (λ, u) and p t 0 (λ, u) for a single generalized ray. The difficulties of rebinning of the offset multi-source geometry rays to a single-sheet surface are illustrated in figure 5(c). Due to the positive angle that such rays make with the transaxial plane, we are only able to simultaneously fit a half of each the ray (λ, u, V (λ, u)) and its complementary ray (λ c , u c , V (λ c , u c )) to a single-sheet surface (e.g. the bottom sheet) while the other half diverges, unavoidably resulting in a large axial deviation for the singlesheet surface rebinning of offset multi-source geometry data. We propose to circumvent this problem by way of rebinning to a multi-sheet surface. In the following we make use of the mixed 2D fan beam transform on the multi-sheet surface ζ 0 defined as Note, that the first integral is taken over ζ b 0 and the second over ζ t 0 (see also figure 5(b)) thus each generalized ray breaks between the sheets. Which sheets the integrals are taken over is determined by the mapping of ray segments 1 = [l 0 (u) − l(u), l 0 (u)], 2 = [l 0 (u), l 0 (u) + l(u)] to sheets, μ : {1, 2} → {b, t} illustrated in figure 6(a). Here we chose to map the first half of the ray 1 to the bottom sheet, ζ b 0 (μ(1) = b), and the second half 2 to the top sheet ζ t 0 (μ(2) = t). For two-sheet surface, there are only two possible choices of such mapping, both resulting in the same multi-sheet surface, but with the labels of the sheets permuted.
Though in general this mixed fan beam transform cannot be measured exactly using line measurements, as visualized in figures 5(b) and (c) it matches the data measured by an offset multi-source geometry far better than any 2D transform on a single-sheet surface. The mixed fan beam transformg 0 (λ, u) can hence be approximated by g 0 (λ, u), the rescaled axially truncated cone beam data (scaling factor as in (18)) rebinned to the multi-sheet surface ζ 0 using the rebinning function V 0 In the next section we describe how to choose the multi-sheet surface ζ 0 and the rebinning function V 0 such that this approximation is optimal with respect to a chosen quality measure Q(ζ 0 , V 0 ). We now derive the relation between the fan beam transform on a two-sheet surface and the mixed fan beam transform. To this end we split the formula in equation (22) into integrals along halves of the generalized ray. Then (22) u, l)), s = μ(i) and recall 1 = [l 0 (u) − l(u), l 0 (u)], 2 = [l 0 (u), l 0 (u) + l(u)]. With this notation it holds for the fan beam transform on each sheet where (λ c , u c ) denotes the complementary ray to (λ, u). While p s 0 is exactly what we would like to approximate (we can deal with a transform on a single-sheet surface), the right hand side involves integration along half rays. Therefore it cannot be approximated with the measured data and so (24) seems of no immediate use. However, substituting (24) into the definition of the fan beam transform on the multi-sheet surface (21) yields , which after regrouping terms we see is identical to Equation (25) lies at the heart of the MSSR. It relates the fan beam transform on the two-sheet surface to the mixed fan beam transforms corresponding to the primary and complementary rays. In turn, each of the mixed fan beam transforms on the right hand side can be approximated by the offset multi-source geometry data g 0 (λ, u) and g 0 (λ c , u c ) providing an approximation to the unknown fan beam transform on a two-sheet surface p 0 (λ, u) This approximation is graphically depicted in figure 5. Notice, that in contrast to the single-sheet surface rebinning, equations (25) and (26) give an implicit approximate relation between the fan beam transform on the sheets p s 0 and the rebinned data g 0 .

Rebinning equation for the multi-sheet surface
The principle can be generalized to a surface with any even number of sheets S, ζ s 0 , s = 1, . . . , S with the objective to extend the method to trajectories with winding numbers t > 1. An example of such a trajectory is t-threaded helix, where the threads are not equi-angularly distributed (note that such trajectory is still invariant under the simultaneous rotation and translation θ ). In particular for a trajectory with winding number t, the multi-sheet surface constructed will have S = 2t sheets. Let 1 , 2 , . . . , S be a partition of the planar ray (λ, u, l), l ∈ = [l 0 (u) − l(u), l 0 (u) + l(u)] along its l coordinate into closed subintervals i , i = 1, . . . , S : ∪ i i = , ∩ i˚ i = ∅. Furthermore, we are going to restrict the partition to be symmetric with respect to the centre of the ray l 0 (u) i.e. it holds S−i+1 = 2l 0 (u) − i , i = 1, . . . , S/2, where the negative sign denotes the reflection of the interval around 0. This ensures that the primary and complementary rays are treated in the same way, retaining the symmetry of the underlying planar problem. For S = 2 (two-sheet surface) such partitioning corresponds to splitting the ray in halves.
To accommodate rebinning of multiple cone beam rays per generalized ray we extend the mapping μ as follows where μ i (r) := μ(i, r) = s such that in each interval i ⊂ , the ith ray segment (λ, u, v, l), l ∈ i (u) of the r(λ)th cone beam ray (λ, u, v) is assigned to the sheet ζ s 0 . The choice of μ determines which i patches will be grouped together into one sheet. The optimal rebinning surface and function pair also nontrivially depends on the partition of the interval .
Rebinning t cone beam rays to one generalized ray necessitates t mixed 2D fan beam transforms on the surface ζ 0 (one for each cone beam ray) Note, that the sheet on which f is integrated along the ray segment i is given by the mapping μ i (r). The rth mixed fan beam transformg (r) 0 (λ, u) has been constructed so that it can be approximated via optimal rebinning with the offset multi-source geometry data , but whenever convenient we consider it as a multi-valued function V r(λ) Analogously to the two-sheet case, we obtain the relation between the fan beam transform on the multi-sheet surface p 0 , and the t mixed fan beam transforms (29)

Optimal multi-sheet surface and rebinning function
In [4] we considered a constrained optimization problem for a single-sheet surface and rebinning function for the offset multi-source geometry. The constraints on the detector vertical extent, result in the constraint on the rebinning function and they have been accommodated using Lagrange multipliers ν 1 , ν 2 . Here we adapt the same approach for finding the optimal multi-sheet surface and rebinning function. The objective is to find a pair (ζ 0 , V 0 ) which minimizes the following axial deviation functional where where the bounds v 1 , v 2 have been periodically extended from [−π, π ) to [−π, (2t − 1)π ).
For q > 1 the problem is strictly convex, the box constrains are qualified, thus (32)-(34) has a unique solution.
Using Lagrange multipliers we arrive at the following unconstrained cost functional In principle any positive bounded weights w will lead to a problem with unique solution. However, we are going to choose the weights so that they are consistent with the underlying backprojection formula. To discuss the choice of weights w we reparametrize the integral (35) U (λ, x, y), L(λ, x, y) where i is defined by L(λ, x, y) ∈ i and we made use of the Jacobian of the change of coordinates between (u, l) and (x, y), By construction in the volume integral (35) each ray on each sheet is integrated along exactly once (without redundancy). Therefore the weights W (λ, x, y)/L(λ, x, y) are chosen as in nonredundant fan backprojection formula (16) derived in section 3.2 and hence Going back to the coordinates (u, l), we have For a system which is not θ -symmetric we can formulate the same optimization problem independently for each of the surfaces and the corresponding rebinning function (ζ k , V k ), associated with the rebinning centres λ k , k = 1, . . . , N c within one period of the effective trajectory.
Equation (36) considered for a single sheet at the minimum of the cost functional gives rise to a point by point way to evaluate the fit between the rays and the sheet ζ s where we use a shorthand notation i.e. the integration is only carried over ray segments rebinned to the sheet ζ s 0 . An appropriate normalization for the axial deviation map Q s q is the volume The last step follows using the reparametrization of the integral with complementary variables as in section 3.2. Note, that this is exactly the same set of weights and normalization, which is applied to ζ 0 in (45). At each point (x, y), Q s q (x, y) gives the L q axial deviation of the rebinned rays from the sth sheet, ζ s 0 , of the multi-sheet surface. In the second paper in this series we are going to show how we can actively exploit the information of the ray to surface fit to derive models of the axial blur for the deconvolution step discussed in section 8.
In the following we consider q = 2 because in this case the unique minimizer can be found by means of a globally convergent iteration. However, other choices of q are plausible e.g. the choice q = 1 would result in a penalty directly proportional to the ray divergence. In principle, once the solution of the least squares problem can be obtained, any q 1 norm solution can be determined with for example iteratively reweighted least squares method [19].
In the particular case of q = 2 the objective functional is a weighted square axial deviation of the rays from the surface Setting the gradient ∂Q ∂(ζ 0 ,V 0 ) = 0 and solving for ζ 0 and V 0 we obtaiñ and was defined in (41).

proof in the appendix for definition).
Proof. See the appendix.

Examples of multi-sheet rebinning surfaces and functions
In this section we give examples of multi-sheet surfaces and rebinning functions for different trajectories. We restrict ourselves to the trajectories with period 2π , t = 1, and two-sheet surfaces. In our experiments we used an offset multi-source geometry with a static ring containing 1152 sources and 10 rings each of 1400 detectors. The offset of the first detector ring from the plane of sources was 10 mm and the spacing of the detector rings 2 mm. The radii of the source and detector rings are 600 and 450 mm, respectively and the radius of the field-of-view, R FOV , is 400 mm. We assumed a constant axial translation of the scanner of 16 mm per period of the effective trajectory, 2π .

Optimal two-sheet pair
7.1.1. Helical trajectory. We first compute the optimal pair for a helical trajectory. Figure  7 shows the optimal two-sheet surface ζ 0 and the corresponding rebinning function V 0 . Along with the surface, we plot the effective source trajectory subsampled by factor 4 for visualization purposes. V 0 has been binned to the corresponding detector rows on the cylindrical detector. For each source within the 2π range of projections rebinned to ζ 0 , for each of the active detectors, the image V 0 colour codes the corresponding row on the cylindrical detector.
In figure 8 we plotted the normalized weighted L q axial deviation map (40) per pixel, (Q s q (x, y)/C(x, y)) 1/q . The first row corresponds to the L 1 -norm, Q 1 , and the second to the L 2 -norm, Q 2 . The shape of the axial deviation maps is quite asymmetric, which is due to the heavy asymmetry of the helical trajectory. While Q 1 is a more robust and more quantitative measure of the axial deviation, Q 2 shows the finer details of the rays to surface fit.

4-threaded helix trajectory.
In the second experiment we used a 4-threaded helix as a trajectory. The optimal rebinning surface and function are shown in figure 9. The rebinning surface is more regular than the surface corresponding to a single helix. However, for the axial deviation maps in figure 10, the asymmetry is still pronounced. In particular, we observe that the axial deviation seems to become larger around the jumps of the effective trajectory, and the region of relatively low axial deviation is cross shaped following the 4-threads of the helix.

32-threaded helix trajectory.
In the final example, the number of helices is close to √ N s [22], resulting in the setup with the most radial symmetry. Figure 11 shows the optimal rebinning pair. The rebinning surface looks very similar to the double-cone surface, a parametrized version of the MSSR method developed in [2]. The corresponding axial deviation maps Q 1 , Q 2 in figure 12 are almost radially symmetric.

Source firing sequence
The problem of finding an optimal source firing sequence (a sampled effective trajectory), in the sense of minimizing square axial deviation is a combinatorial problem. Hence, we restrict ourselves to comparing the axial deviation values of the here tested trajectories: single helix, 4-threaded helix and 32-threaded helix. Table 1 shows the values of Q 1 , Q 2 corresponding to the L 2 optimal rebinning pairs for those trajectories (recall that the rebinning pairs for each trajectory were optimized w.r.t. Q 2 ). While the single helix has the smallest axial deviation (in both L 1 and L 2 norms, Q 1 , Q 2 ), and the axial deviation seems to increase with the number of threads, this increase is not significant, and is outweighed by the benefits of the close to radial symmetry of the rebinning surface at the deconvolution stage [2].

Axial deconvolution
In section 5.3 we obtained an implicit relation between the fan beam transforms on the sheets of the surface, p s 0 , and the data rebinned to multi-sheet surface g 0 , (25), (26). To apply standard rebinning method one would first attempt to separate the transforms on each of the sheets p s 0 , reconstruct each of them individually and finally interpolate the resulting volumetric image from the images on the individual sheets (these are simple surfaces) that intersect the bounded cylindrical region containing the support of f . Unfortunately (21) is severely underdetermined and hence it is not possible to solve it without additional information. Furthermore, as the data on any two different multi-sheet surfaces is not correlated in a simple way, there is no plausible way to set up a meaningful deconvolution problem in the data space.
However, the linearity of the Radon transform and hence of its inverse (filtered backprojection, B) enables us to shift the deconvolution of the information on the individual sheets of the multi-sheet surface from the data to the image domain, where there is an evident correlation. For a multi-sheet surface ζ m we have where M Z := {λ m : ζ λ m ∩ = ∅} is the set of rebinning centres of all the multi-sheet surfaces intersecting the bounded cylindrical volume .
As the multi-sheet surfaces filling the volume intersect, the equations (46) are not independent but constitute a system of coupled linear equations defining the axial deconvolution problem for unmixing the information from multi-sheet surfaces. The system (46) can be uniquely solved provided the multi-sheet surfaces sample the volume densely enough. An important feature of the system (46) is that it decouples into independent blocks for each fixed value of (x, y), exactly as in the case of axial interpolation employed in single-sheet surface rebinning methods. This and other properties of the deconvolution problem (46) allow us to devise efficient numerical procedures for its solution discussed in the second paper in this series [3].

Conclusions
We introduced the idea of rebinning to multi-sheet surfaces, which is a basis for a new family of rebinning methods for reconstruction from axially asymmetrically truncated 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 the first paper of this series we have discussed the underlying fan beam transforms and their approximation with the truncated projection data obtained from offset multi-source geometries. We set up a variational problem for optimal rebinning surface and function pair and gave a globally convergent iteration for its solution in the case of a quadratic fidelity term. We gave examples of different optimal rebinning pairs for different trajectory choices with period 2π . Due to rebinning to multi-sheet surfaces, the reconstructed 2D problems define the volumetric image implicitly. We concluded formulating an axial deconvolution problem for its recovery. In the second paper in this series we discuss details of the 2D reconstruction on multi-sheet surfaces including strategies for dealing with approximated data. We introduce improved deconvolution models using distribution of the axial distance of the rays to the surface and explain how the deconvolution step can be efficiently and robustly implemented. Finally, we demonstrate the performance of the methods on both the simulated data and the real data collected with Rapiscan RTT80 scanner. grateful to Dr Edward Morton (who invented the RTT system), Ken Mann and the rest of the RTT80 development team at Rapiscan who have given us access to technical details of the system and data. We would also like to thank the anonymous reviewers for their invaluable suggestions for improvement of the structure of both parts of the manuscript.

Appendix. Proof of theorem 1 (convergence of the alternating iteration)
We show that the alternating iteration (44), (45) converges to the unique minimum of the Lagrange function (43). The proof is motivated by the proof in [4] for the constrained surface problem, which itself is based on the proof for the unconstrained surface problem in [8]. In the sequel we need the definitions of the following function spaces: • A space of admissible images on a multi-sheet surface: a weighted Lebesgue space is the L 2 W ( 2D )-norm we use on an individual sheet f s and the weight where is the integral taken over λ : L(λ, x, y) L 0 . (A.2) is a non-redundant fan beam backprojection of the unit data and it is positive everywhere in the interior of 2D . We used a shorthand notation L = L(λ, x, y) and U = U (λ, x, y).
• A space of 2tπ fan beam transforms: a standard Lebesgue space L 2 (Z 2D,t ) where Z 2D,t := [−π, (2t − 1)π ) × [−u max , u max ] with the induced norm We define a linear operator Up to a multiplicative factor, P is a weighted Radon transform of a multi-sheet function f . Hence, P is composition of a bounded operator (weighting with bounded weights 0 < 2l 0 /(l 0 + l) 2l 0 /(2l 0 − l) 2l 0 /(l 0 − l) < ∞) and a compact operator (Radon transform) and as such it is a compact operator. Hence for a bounded domain, P has a discrete spectrum. Using the formal definition (g, P f ) = (P * g, f ) W , we obtain the adjoint operator P * P * : L 2 (Z 2D,t ) → L 2 W (S × 2D ) g(λ, u) → (P * g)(x, y) Recall that (41) In the unconstrained case, as considered in [8], which in turn for a compact operator P implies that the largest eigenvalue of P * P, σ 2 1 < 1 and hence P * P W < 1.
Proof. As in the mixed fan beam transform, each ray is represented exactly once, we cannot use dual parametrization of the same line (the primary and complementary rays), which would easily provide the argument. Instead, we are going to rely on the fact that each point P has infinitely many representations in the ray coordinates, i.e. for each angleλ there exist the different corresponding values of u and l.
Let O, P, P be points which belong to the same sheet f s O : (λ a , u a , l a ), (λ b , u b , l b ) P : (λ a , u a , l a,P ), (λ P , u P , l P ) P : (λ b , u b , l b,P ), (λ P , u P , l P ) with l a , l b , l P , l P , l a,P , l b,P ∈ i and r(λ a ) = r(λ b ) = r(λ P ), see also figure A1. Let f s ( (λ a , u a , l a ), R(λ a , u a , l a )) = c, where (λ a , u a , l a ), R(λ a , u a , l a )) are the polar coordinates of O and c is some real constant.