Tomography with extended sources : Theory , error estimates , and a reconstruction algorithm

Recently, an approach to tomography with extended anisotropic radiation sources has been introduced, which helps to overcome the challenges resulting from the low brilliance typical for x-ray laboratory sources. The method is based on the three-dimensional Radon transform (3DRT) which uses planar integrals instead of line integrals. By extending the source spot in one direction, more photons can contribute to image formation while the impact on the resolution is minor with the 3DRT approach. In this work we present a more comprehensive description of the method, derive quantitative error estimates for the extraction of these planar integrals measured with a finite source size, and validate the 3DRT scheme by analytical theory. We also demonstrate a simple and efficient reconstruction algorithm for 3D Radon data. Finally, we further substantiate the method with experimental results obtained at a microfocus x-ray source with an extremely anisotropic source spot.


I. INTRODUCTION
X-ray computed tomography (CT) is widely used today for nondestructive three-dimensional (3D) imaging in medical [1], biological, and material science applications [2].Apart from the high penetration power, hard x rays also offer further advantages: The broadening of the beam by scattering and diffraction within the object are sufficiently low so that the projection approximation required for CT holds up to large specimen thickness and small voxel size.Moreover, the absorption coefficient is quantitatively tractable and can be linked to material composition.Finally, the small wavelength holds the promise of reaching high resolution.In practice, however, most applications of analytical CT reach at best moderate resolution in the micron range, while nanoscale resolution below that of visible light remains rare.This is primarily due to the limited brilliance of laboratory x-ray sources, despite some important advances such as the liquid jet anode [3][4][5][6] or the inverse Compton scattering based sources [7,8].For this reason, synchrotron radiation (SR) is still indispensable for tomography at the nanoscale.X-ray phase contrast CT based on free space propagation [9][10][11], which enables high resolution also for non-or weakly absorbing soft tissues, are even more constrained to SR.In order to reach the required degree of spatial coherence or the resolution, one cannot simply slit down or refocus a laboratory x-ray source, without detrimental loss of flux.
To overcome these constraints, we have recently proposed an approach to analytical CT based on the 3D Radon transform (3DRT) [12].The 3D Radon transform is well understood theoretically [13] and several reconstruction methods for conebeam data based on the 3DRT have been considered (see, i.e., Refs.[14,15]).The method put forward in Ref. [12] is designed to make a larger fraction of photons exploitable without losing resolution or coherence by relaxing the brilliance or source size Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
in one of the two source dimensions.Extending the source in one direction s x while keeping the other s y sufficiently small (as required for resolution or coherence), a larger flux can contribute to the image.At the same time the "good" image quality (resolution, coherence) as determined by the small source direction s y is translated isotropically into all (3D) directions of the reconstruction volume.This 3DRT approach requires the measurement of planar integrals, i.e., integrals over planes through the object, instead of using line integrals of the conventional 2D Radon transform.Further, the normal vector of the planes have to be sampled on the unit sphere.Due to constraints imposed by the anisotropic source, two rotation axes have to be used.While we have previously shown by both simulation and experiment that the 3DRT method works in principle [12], we had not yet given quantitative limits on the tolerable source extension s x s y with associated error bounds, which is the scope of the present work.
In this work, we hence consider to which extent detector line integrals, measured in cone-beam geometry with a source of given width s x , approximate the planar integrals required for the 3DRT.Note that the error studied here by analytic theory includes but is not limited to the well-known cone-beam error controlled by the cone-beam angle σ , and depends both on σ and the source size s x .The results are useful also for the conventional case of micro-CT with isotropic source size.As we discuss in detail, two conditions have to be distinguished: First, the nonlinear averaging of the signal such as in absorption contrast, where the negative logarithm of the detector value has to be taken prior to reconstruction.Second, the linear averaging which applies, for example, for phase contrast or x-ray fluorescence tomography, but also as the limit of weakly absorbing samples.In the linear case, we show that the finite source size results only in a minor error contribution with no loss of resolution.Contrarily, for the nonlinear averaging associated with absorption contrast, the tolerable extension of the source in s x is more critical and must be chosen according to object properties-in particular the spatial variation of the projected absorption coefficient.
We then show experimentally that sharp reconstructions of a biological object (gerbil cochlea) can indeed be obtained based on the 3DRT, resulting in isotropic resolution on the order of 20 μm, even when projections are almost completely blurred FIG. 1. Sketch of the cone-beam geometry.A point source at r 0 projects some density of the object, described by f , on a twodimensional pixel detector.A line L on the detector defines a plane (red) in R 3 together with the source.Integration along a set of parallel lines yields a one-dimensional profile which approximates the planar integrals.Each line, defining a plane, corresponds to one value of the profile.A full set of planar integrals can be used to reconstruct the original density f by the inverse 3DRT.
in one direction owing to a highly anisotropic source size of 0.04 mm × 12 mm.With respect to the previous experiments [12], the angular sampling, data processing routines, and the reconstruction procedure were improved further.
The manuscript is organized as follows: Following these motivational paragraphs, the cone-beam reconstruction scheme and basic definitions are introduced.Section II presents the analytic theory, including the geometric errors of the linear and nonlinear case, as well the related quantification of the resolution loss.An explicit reconstruction algorithm is given in Sec.III, addressing also the issue of angular sampling for the 3DRT.Finally, the experimental demonstration is presented in Sec.IV, before the manuscript closes with a brief conclusion.

A. Cone-beam reconstruction scheme
The cone-beam geometry (Fig. 1) is of outstanding practical importance for computed tomography: it naturally arises, whenever projections from a point source are measured with a two-dimensional camera.These cone-beam projections can be defined as where r 0 and r d denote the position of the source and a point on the detector, respectively.The function f describes some local property of the sample, for example, the attenuation of x rays of a specific energy per unit length.A single acquisition naturally produces p(r 0 ,r d ) for all r d on the detector.The original function f can be reconstructed exactly from a set of cone-beam projections from different directions, if the set of directions satisfies certain completeness conditions [16,17].
Varying the direction can be achieved by either rotating the sample and thus f , or by moving a gantry with source and detector on an orbit around the sample.
Three-dimensional imaging in cone-beam geometry has been studied extensively and several reconstruction algorithms have been proposed [1,[15][16][17][18].One approach considered in this work is based on the three-dimensional Radon transform (3DRT), ( It maps a function, f : R 3 → R + , to the set of its planar integrals, g : S 2 × R → R + .An explicit inversion formula for the Radon transform is given by [13] f Here "•" denotes the usual inner product on R 3 , and n ∈ S 2 is understood as a unit vector in R 3 .The integration is computed over the unit 2-sphere S 2 with integration variable n and dS denotes the surface element on the 2-sphere.To relate the conebeam projections to the Radon transform recall the following basic geometric fact: In cone-beam geometry, any line L on the detector together with the position of the source r 0 determines a plane in R 3 .As a direct consequence, planar integrals can be extracted up to a non-negative kernel K pnt by integrating p(r 0 ,r d ) along a line L on the detector, We suppress the relation between the Radon coordinates (n,s) and r 0 and L parametrizing the same plane for now.The purpose of the weighting function w L is to compensate the missing curvature of the detector.The explicit form of the kernel K pnt and the choice of the weighting function w L will be discussed in Sec.II C. The cone angle σ , the angular field of view of the object seen from the source, constitutes an important parameter of the setup.It can be defined as σ := max{ (r − r 0 ,r − r 0 ) : r,r ∈ supp f }.
In particular, we have that K pnt → 1 for σ → 0. Hence, for σ 1, the approximate planar integrals defined in Eq. ( 4) can be used to reconstruct f .This provides us with a particularly simple method for cone-beam tomography with small cone angles: (1) Measure cone-beam projections p for various orientations.
(2) Extract (approximate) planar integrals g by integrating along (straight) lines on the detector.
(3) Reconstruct f from the planar integrals using some inversion algorithm for the three-dimensional Radon transform.
In this work we analyze the compatibility of this approach with an extended source.

A. Optical system
Consider the geometry as sketched in Fig. 2. We describe the object's local attenuation coefficient by f (r).Let the object be contained in a ball of radius R around r 1 , FIG. 2. Cone-beam geometry with an extended source.The object, described by f , is contained in a ball of radius R around r 1 .The source is centered at r 0 .The detector lies in the plane through r 2 with normal z (r 2 − r 0 ).The source to object distance r 1 − r 0 and the object to detector distance r 2 − r 1 are denoted by z 01 and z 12 , respectively.If the (line) source S is parallel to a line on the detector L, all rays from points r s ∈ S to points r d ∈ L lie in a common plane P (light red).
supp f ⊂ Ball(r 1 ,R).Consider a fully incoherent and isotropically emitting source.The latter assumption is justified if we restrict our consideration to a narrow cone around the optical axis z.We describe the source by its (spectrally and angular integrated) radiant emittance I 0 j (r s ) in the plane through r 0 with normal z (z = z 0 ).For simplicity we assume j to be non-negative and normalized such that j (r s )dr s = 1.Naturally, the source spot is assumed to be finite with spatial extent s x and s y in x and y directions, respectively, thus supp j = S ⊂ {(x,y,z 0 ) : |x| < s x ,|y| < s y }.
By neglecting any diffraction within the object we can write the total intensity in the detection plane as where p was defined in Eq. ( 1).In the empty beam without an object we have p ≡ 0 and the intensity is For relatively small sources we have Division by the empty-beam intensity then yields For a point source, the negative logarithm of u resembles p exactly.In the general case we obtain the effective projection, With a finite source spot we have that p = p and the image is corrupted.

B. Extraction of planar integrals from blurred projections
For our further analysis we assume that s y s x , such that x is the low-coherence direction and y the high-coherence direction.For now, we consider the limit s y → 0, i.e., let S be a line.As we will see, s x > 0 does not affect the resolution of our method but rather introduces a systematic error.The general case, s y > 0, will be considered in Sec.II E.
Let L be a line on the detector, the integration line, and define where M = z 02 /z 01 is the geometric magnification and w L some weighting function.Inserting Eq. ( 6) into Eq.( 7), we obtain Note that Eq. ( 8) is the generalization of Eq. ( 4) for a line source.Again, g only depends on the values of f in the set of planes defined by L and the source points r s ∈ S.More precisely, g is determined by f | conv(S∪L) where "conv" denotes the convex hull in R 3 .Accordingly, if S L (as lines in R 3 ) then S and L lie in a common plane P with conv(S ∪ L) ⊂ P and g is determined by f | P .This simple fact allows us to compare g to the corresponding planar integral, In fact we could write g = P[f | P ] and g = P[f | P ] and study the operators P,P : L 2 (R 2 ) → R. The constraint S L ensures that P also operates on L 2 (R 2 ) instead of L 2 (R 3 ) and thereby makes P and P comparable.However, we are not going to use this notation.For the following analysis, we require L S and hence, for our choice of coordinates, L x.

C. Geometric errors
We first consider the limiting linear case, where the exponential and logarithm in Eq. ( 8) can be linearized and This linearization can be made for a weakly absorbing object where p 1, or more generally, for any process where linear projections are measured.Substituting Eq. ( 1) into Eq.( 9) and interchanging the order of integration we obtain Here we introduce x s and x d as the x ordinates of r s = (x s ,0,z 0 ) ∈ R 3 and r d = (x d ,y d ,z 2 ) ∈ R 3 , respectively.For fixed r s , the integration over dt and dx d can be transformed to a surface integral dA over a subset of P .In the following we are going to derive an explicit expression for this surface integral.This will ultimately enable us to compare Eq. ( 9) to the planar integral g.As a first step, we define the auxiliary polar coordinate θ as For its differential we obtain the relation, with r s fixed.Note that (θ,t) constitute polar coordinates around r s in the plane P .The area differential in the plane consequently reads In order to make the right-hand side of Eq. ( 11) independent of r s and r, we employ the approximations r − r s ≈ z 01 and Note that r d on the left-hand side of Eq. ( 12) is not an independent variable but the linear projection of the point r from r s onto the detector plane.The natural choice of the weighting function can be read of from Eq. ( 12), From here on we will suppress the subscript L in w for notational simplicity.Under the assumption that the blurred projection is fully captured by the detector, we can formally extend the integration to the entire plane P .Inserting Eqs. ( 12) and ( 13) into Eq.( 10) and exchanging the order of integration a second time yields where we introduced the integration kernels, Equation (14a) is independent of the source size and purely caused by the diverging beam in cone-beam geometry.In addition, the finite source causes the second kernel, Eq. (14b), which reduces to 1 for a point source.By expanding the kernels into Taylor series in r around r 1 , we find for r ∈ Ball(r 1 ,R), where m k is the kth moment of the emittance distribution j , m k = S j (r s ) r s − r 0 k dr s and σ = R/z 01 the cone angle.This leads to the upper error bound, The moments of the emittance distribution are limited by the source dimension, m k s k x , hence m k /z k 01 (s x /z 01 ) k .The first moment m 1 can be eliminated by adapting the coordinate system to the center of mass of the source distribution j .For most applications, one can certainly assume s x /z 01 1 in which case the finite source is negligible and the right-hand side of Eq. ( 15) is effectively given by σ .

D. Nonlinear averaging
We now consider the nonlinear case in order to study the nonlinear averaging which takes place in Eqs. ( 6) and ( 7).This effect is in itself independent of the cone-beam geometry.In order to analyze it independently we make the following assumptions: (i) small cone angle σ ≈ 0 and hence w = 1; (ii) a thin object such that the projection integrals p(r s ,r d ) can be approximated by functions h : R 2 → R, with Assumption (ii) is equivalent to approximating the threedimensional density f with a two-dimensional screen at z 1 .
The argument of h on the right-hand side of Eq. ( 16) gives the intersection of rays from r s to r d with this screen.Inserting Eq. ( 16) into Eq.( 5) yields where Here * denotes convolution with respect to the argument marked with "•".This j eff describes the emittance rescaled to the object plane z = z 1 .The effective projection can be expressed as where we defined the exponential convolution * logexp .From now on, we will suppress the second argument y d /M of h.
Using this notation, we may write the approximate planar integral as with L := M −1 L. The difference to the exact planar integral g then reads We denote by (x d ) the integrand of Eq. ( 18), Equation ( 19) quantifies the local approximation error.Writing the two convolutions explicitly and pulling the second term into the logarithm yields where ) is the deviation from the local average of h in the vicinity of x d .In the second line of Eq. ( 20) we used that multiplication with a constant commutes with convolution.Equation (20) gives us a lower bound for , 0 (x d ) − ln cosh (δh(x d )), where δh(x d ) = max x∈S | x d (x)| quantifies the maximal variation of h within the effective source size.The effect of the exponential convolution is shown in Fig. 3. There, a profile h featuring sharp edges is plotted together with j * logexp h, and j * h where j is the emittance distribution of a uniform source, Eq. ( 23).In the lower half of Fig. 3 the relative deviation a = (j * logexp ah − j * ah)/a is plotted for different a ∈ R + .It is strictly negative and its modulus increases with increasing a.Note that for 5% transmission, the relative deviation reaches 20% in this example.
If | x d | is sufficiently small, we can expand Eq. ( 20) into a Taylor series.The first nonzero term is By integration of Eq. ( 21), we obtain Equation (22) gives the lowest order of the error due to the exponential convolution.The assumption of weakly varying absorption (| x d | 1) which led to Eq. ( 22) breaks down for The numbers a = 1, 3, and 5 correspond to approximately 37%, 5%, and 0.7% transmission, respectively.sharp edges in the projected attenuation.We analyze this case by considering a piece-wise constant function, We assume that the effective source size s x = s x M/(M − 1) is smaller than the structure size of the object, such that a i+1 − a i > 2s x for all i.For simplicity, we assume the source to be uniform, By carrying out the integration in Eq. ( 17), we get , As a numerical example we consider a rectangular function (N = 1) of width a 1 − a 0 and height h 0 , giving The second fraction on the right-hand side of Eq. ( 24) is asymptotically linear near 0 and bounded by 1.Note that h 0 = − ln μ where μ is the projected attenuation coefficient.For ≈ 5% transmission, we have h 0 ≈ 3 and (h 0 coth h 0 − 1)/h 0 ≈ 0.67.In this case |g logexp − g|/g ≈ 1.34(2s x )/(a 1 − a 0 ).If s x is of comparable magnitude as the structure size a 1 − a 0 , the relative error is of the order 1.This exponential error appears in conventional CT as well.In medical imaging it is one source of metal streak artefacts and also known as exponential edge-gradient effect (EEGE) in the literature [19].It is, however, not as apparent there, as the source size limits the resolution and is thus usually kept smaller.For our method, without any other restriction on the source size, this nonlinear error does impose a limit on s x .

E. Resolution
We have seen that the one-dimensional (line) source does not influence the resolution if the integration line is parallel to the source.Instead, it affects the quantitativity of the planar integrals by a geometric and a nonlinear contribution.We gave upper bounds (up to higher orders of the cone angle) for these errors and showed that the geometric contribution is negligible whereas the nonlinear contribution is determined by the strength of local variations of the projections.
Let us now consider the general case where 0 < s y s x .Let r s = (x s ,y s ,z 0 ).Note that the emittance distribution j is a two-dimensional density here instead of a one-dimensional one.First, we analyze the linear case given by Eq. ( 9).We decompose S into the set of parallel lines with constant y s , defining S y s := {(x s ,y s ,z 0 ) : x s ∈ R} ∩ supp j .Further, define the one-dimensional normalized densities j (y s )dy s := S ys j (x s ,y s )dx s dy s , j y s (x s )dx s := j (x s ,y s )dx s dy s j (y s )dy s .
For fixed integration line and sample orientation, each of these lines S y s on the source satisfies S y s L and for this reason fixes a (different) plane.These planes can be parametrized by their Radon coordinates (n(y s ),s(y s )).Following this approach, the generalization of Eq. ( 9) for a two-dimensional source reads where FIG. 4. Rectangular support for different orientations of the coordinate system.The numbers s x and s y are defined as the smallest number such that the rectangle [−s x ,s x ] × y ,s y ] (with respect to the coordinate system defined by x) fully contains S (red).If x is chosen such that s y is minimal, any other direction x results in a bigger sy .If S is unknown, only an upper limit to sy as a function of (x,x), s x , and s y can be computed.
The planar integrals are blurred with the point spread function (PSF) j and Eq. ( 27) gives the extent of the blurring.Importantly, the blurring depends only on s y , and not on on s x .For the nonlinear case, Eq. ( 8), we write g logexp = g w + (g logexp − g w ).The term in parentheses is given by Eq. ( 18), the convolutions being two-dimensional in this case.The nonlinear generalization of Eq. ( 26) then reads with (x d ) as defined in Eq. ( 19).The exponential convolution further distorts the measurement in addition to the blurring.The properties of this nonlinear contribution were discussed in Sec.II D.
It may not always be experimentally possible or desirable to exactly align the integration line L with the low-coherence direction.Some source shapes do not even have a unique lowcoherence direction.We may define s x and s y to be the smallest real numbers such that supp j is contained in the rectangle with sides 2s x and 2s y which are aligned with the coordinate axes.Let x be a direction such that s y is minimal and let x be any other direction that is rotated by the angle α with respect to x.The source dimensions in this rotated frame are bounded by (cf.Fig. 4) For α 1 this gives Since s y effectively limits the resolution by Eq. ( 27), this measures how fast the resolution deteriorates, if the integration line is not strictly aligned with the low-coherence direction.
We see that as long as αs x /s y 1 the resolution does not deteriorate significantly.

A. Reconstruction algorithm
For n ∈ S 2 ,s ∈ R let g(n,s) be the values of the Radon transform R of f , as defined in Eq. ( 2).The original function f can be recovered exactly from the full Radon transform Rf , if it fulfills certain smoothness properties [20].
Let ⊂ S 2 + be a discrete subset of the hemisphere with # = N, ⊂ R with # = K, and let G = {g(n,s) | n ∈ ,s ∈ }.We state the discrete reconstruction problem as follows: Given G, approximate f .The number K is effectively given by the number of detector pixels in the high-coherence direction and the number N depends on the chosen sampling scheme.
For notational clarity, we motivate the algorithm starting from the inversion formula Eq. ( 3).See Refs.[1,20] for a derivation of the inversion formula in n dimensions.First, we rewrite Eq. ( 3) to where ĝ is the filtered Radon data, By exploiting the radial symmetry of g, g(n,s) = g(−n, − s), which is conserved by the filtering, the integration domain can be reduced to the hemisphere S 2 + , To discretize the integral in Eq. ( 28), we replace ĝ(n,r • n) with the nearest point in , ĝ(n,r • n) ≈ ĝj (r • n j ) where n j = argmin n j ∈ dist(n,n j ) and ĝj : R → R + , ĝj (s) = ĝ(n j ,s).
Here "dist" denotes the distance on the sphere.This yields a semidiscrete version of Eq. ( 28), where a j = A j 1dS and A j = {n ∈ S 2 + | dist(n,n j ) dist(n,n j )∀j = j }.Note that the A j represent the partitioning of the hemisphere into the Voronoi regions (cf.Ref. [21]) from the generating set .If the n j are distributed uniformly on the hemisphere this simplifies to a j = 2π/N.For an arbitrary sampling scheme the a j can be computed numerically.
In order to use Eq.(29) for numerical reconstruction, the ĝj (s) need to be discretized for s ∈ .One way to do so is nearest-neighbor interpolation ĝj (s) ≈ ĝji = ĝ(n j ,s i ), where s i = argmin s i ∈ |s − s i |.Note that the interpolation here is fundamentally different from nearest-neighbor interpolation in conventional tomographic reconstruction.The method we describe here will instead result in trilinear interpolation.
Here, for each j in Eq. ( 29), the set of s i ∈ corresponds to a set of parallel slabs Here Vol 3 (B ij ∩ V k ) denotes the volume in R 3 .Computation of Eq. ( 30) can be performed for the whole reconstruction volume in parallel.The algorithm consists of two computational steps: (1) For each j , compute ĝj from g j .The filtering can be performed efficiently in Fourier space, using a properly FIG. 5.The Radon coordinates (n,s) are defined with respect to the center of the object r 1 .A line on the detector can be parametrized by its orthogonal direction n 0 in the detector plane and its minimal distance s 0 to the optical axis.The intersection of the optical axis with the detector is denoted by r 2 .Both the Radon coordinates (n,s) as well as the line corresponding to (n 0 ,s 0 ), together with the position of the source r 0 , uniquely describe a plane in R 3 .Their relation can be understood by basic planar geometry.discretized version of F[∂ 2 s g j ](q) = −4π 2 q 2 F[g j ](q), where q denotes the nonangular frequency coordinate.
(2) Back-project the filtered stacks ĝj into the reconstruction volume.

B. Acquisition and sampling
To acquire a full set of Radon data, the sample needs to be rotated with respect to the optical axis (determined by source and detector).We assume that source and detector are mounted on a rigid gantry and describe the rotation by a rotated Cartesian coordinate system (x ,y ,z ) which is related to the fixed sample coordinate system (x,y,z) by θ ∈ SO(3).The detector is parallel to the plane spanned by x ,y .Any integration line L ⊂ span{x ,y } on the detector can thus be parametrized by n 0 ∈ span{x ,y } ∩ S 2 and s 0 ∈ R + such that L(n 0 ,s 0 ) = {r ∈ span{x ,y } | n 0 • r = s 0 }.The plane corresponding to the integration line can be expressed in Radon coordinates (cf.Fig. 5) as where φ = arctan(s 0 /z 02 ) and R(n 0 × z ,φ) denotes a rotation around n 0 × z of angle φ.For small cone angles σ these relations can be approximated by Without constraints on L, all n 0 ∈ span{x ,y } ∩ S 2 parametrize possible integration lines.Every orientation θ consequently yields Radon data g(n,s) for a two-dimensional surface, which for small cone angles becomes a plane.The full three-dimensional Radon space S 2 + × R can thus be acquired by rotation on a one-dimensional orbit 1 ⊂ SO(3), dim 1 = 1.The number N of sampled orientations required for a certain angular resolution θ for this reason scales with The extended source imposes the constraint L x or equivalently n 0 y . (32) Under this constraint every orientation yields Radon data on a one-dimensional curve only, which becomes a line for σ 1.For the full Radon space, rotation on a two-dimensional subset 2 ⊂ SO(3), dim 2 = 2 is required.In this case the number of required orientations N ∼ O(1/ θ 2 ) is squared compared to the unconstrained case.This can be significantly improved by replacing Eq. ( 32) with the relaxed constraint, where α is the maximal admissible rotation angle as discussed in Sec.II E. With Eq. ( 33), every orientation θ yields a strip of width α and rotation on a one-dimensional subset suffices to sample the full Radon space.The number of required orientations is reduced to IV. EXPERIMENT

A. Setup
Experiments were carried out with a sealed tube source (FK-61 04x12 MO, Thomson Tubes Electroniques, France) with a molybdenum anode.The source was mounted with a 6 • tilt to produce a 0.04 mm × 12 mm source spot.A sCMOS pixel detector fiber coupled to a Gadox scintillator (2k × 2k X-Ray SCMOS Camera, Photonic Science, UK) with a resolution of 2048 × 2048 and a pixel size of 6.5 μm was used.The sample was mounted on a Kappa goniometer (515 series, Huber, Germany).A dried cochlea of a gerbil (diameter ≈ 0.5 mm), glued to a pipette tip (Eppendorf, Germany) with epoxy adhesive, was chosen as an object.Note that the object, which was available from a biomedical collaboration [5] does not matter for the scope of this work, but in terms of structure size, composition and absorption is very suitable for the setup.The geometric parameters of the setup are listed in Table I.
Figure 6(a) shows a photograph of the setup.The source was powered by a standard analytical x-ray generator (ID 3003, Seifert, Germany) at 60 kV acceleration voltage with 40 mA anode current.The radiation was pre-hardened by a 180 μm aluminium foil.In total N = 16 001 projections with 0.5 s acquisition time were taken.The different orientations were evenly distributed on the hemisphere by a heuristic algorithm.The total integration time for all projections was 134 min.Each image was dark-field and empty-beam corrected.Nine sets of parallel planar integrals with a maximal angular deviation from  the high-coherence direction of α = 0.5 • were extracted from each projection.Prior to the integration, the top and bottom of each image with respect to the direction of the sample holder were masked out with a Gauß error function to minimize boundary effects.The sets of planar integrals were binned with a factor of 2 in the radial coordinate s.The binned effective pixel size 2p eff approximately matches the effective width of the PSF 2s y .The angular deviation is expected to cause a blurring of less than 32 μm (3.6 times the effective source size) due to the strong anisotropy (s x /s y = 300).Equation (30) was implemented in MATLAB to reconstruct the 3D volume.The set of planes corresponding to one projection was treated as parallel by using Eqs.(31a) and (31b), neglecting the divergence of the cone.The voronoi regions for the discrete spherical integration were computed using the "SphericalVoronoi" algorithm (cf.Reddy et al. [22]) of SCIPY 0.19.1 [23].Reconstruction was carried out on a machine with 2 Intel Xeon E5-2650 v3 CPUs (10 cores each) and 512 GB memory.The computation used all cores and less than 30 GB of the available memory for a 900 3 voxels reconstruction volume and took about 120 min.

B. Results
Figure 6 summarizes the results.The "raw" projection [Fig.6(b)] prior to the reconstruction clearly exhibits a horizontal blurring due to the extended source.The red bar at the top illustrates the effective horizontal extent of the source (in the object plane).Figure 6(c) shows a numerical re-projection of the reconstructed volume for the same orientation.The structure of the sample is sharply resolved in all directions.A slice of the reconstruction volume is shown in Fig. 6(d).It appears to be sharp, as well.Small bone structures are well resolved as can also be seen in the magnified region in Fig. 6(e).Dark cloudy artefacts are visible close to dense bone structures in Fig. 6(d).These may be caused by the nonlinear averaging.After all, these regions absorb up to 75% and thus represent strong edges with variations of h up to 1.4.By measuring the width of different sharp edges with high contrast, we obtained a rough estimate for the resolution of 2-4 px.This corresponds to 20 to 40 μm resolution.It is in good agreement with our estimated blurring due to s y > 0 and α > 0 of 32 μm.
The edges of the glue and the pipette tip show clear signs of edge enhancement due to phase contrast.The bright regions next to the bone edges in Fig. 6(d) may also indicate this.Edge enhancement under similar conditions, although with a much smaller source aspect ratio s x /s y , has already been observed by Vassholz et al. [12].

V. CONCLUSION AND OUTLOOK
We briefly summarize the main error bounds derived for the measuring planar integrals in cone-beam geometry.For the linear case applicable, for example, to phase contrast, the contribution to the total error is bounded according to Eq. ( 15) by O(σ )(m 1 /z 01 ) 2 , where m 1 was the first moment of the source distribution in the "extended" direction, i.e., along x.For analytical CT, we always deal with a small cone angle σ 1 as well as with m 1 /z 01 1.Consequently, this error contribution can be safely neglected for analytical CT.Note that the error remaining for a point source is linear in σ and can be compensated by Grangeat's trick [15].As a result, linear signal integration is compatible with extremely anisotropic sources, for example, s x /s y 10 3 under realistic parameters.Further, we have shown that the resolution in the linear case is not affected by s x , but only the "good" direction s y , since the blurring associated with the source is governed by the length scale s y = s y (M − 1)M −1 [Eq.( 27)], just as in the isotropic case.
Contrarily, for the nonlinear case (absorption contrast) associated with taking the negative logarithm of the detector intensity for each pixel, the validity of approximating the planar integral by the detector line integral can be rapidly compromised when increasing s x .The error, in fact, depends on the variation of the absorption coefficient within the object, more precisely on the variation of the projection integral on a length scale given by the effective source size.For weakly varying absorption, the corresponding error bound is described by Eq. (22).Somewhat more explicit results can be obtained for piece-wise constant projection profiles with arbitrary absorption (or equivalently transmission values); see Eq. (24).Note, however, that the examples given correspond to the quite unlikely case of strong jumps in the profile.Many relevant samples, in particular biological tissues, will exhibit rather smooth transmission.This is because the discontinuities arising from internal boundaries are minimized by integration along the rays through the object, except for very large and flat interfaces.The problem of strong interfaces is well known from artefacts denoted as metal streak artefacts also in conventional CT, although other effects are denoted as metal streak artefacts, as well.More precisely, the described effect of nonlinear averaging is often denoted as the exponential edgegradient effect (EEGE).Indeed, the results on the nonlinear averaging derived here, in particular Eq. ( 20), may also be applied to conventional CT, helping to control the errors and to optimize parameters in instrumental design to avoid such artefacts.
To put the error of the nonlinear case into a practical perspective, in particular for biological tissue, we have carried out an experiment in absorption contrast.It shows that the typical variation of absorption in an object composed of both bone and soft tissue does not impede successful reconstruction with 3DRT, even at extremely anisotropic source sizes.To this end, consider the experimental results illustrated in Fig. 6.Despite a seemingly hopeless blurring of the projections [Fig.6 As we have demonstrated, the herein described 3DRT scheme is compatible with extended sources.As such, it may circumvent the compromise between resolution and acquisition time that one usually has to make for sources of low brilliance.The Appendix contains a simulated comparison between the 3DRT method and conventional tomography based on the 2D Radon transform (2DRT) with filtered backprojection (FBP) for specific parameters.In the simulated case, the 3DRT scheme increases the signal-to-noise ratio by a factor of 4 while the total acquisition time and the surface power density of the source remain equal.
In conclusion, the analytical results derived here and the experimental demonstration of the 3DRT reconstruction in an extremely anisotropic case with s x /s y 300 further validate the method proposed in the recent letter [12].They prove that the 3DRT scheme makes reconstruction compatible with a source size extended along one dimension (the "lowresolution direction"), while the beam properties along the "high-resolution direction" determine the resolution of the entire 3D reconstruction.
Importantly, the current work provides quantitative error bounds useful for experimental design, as well as data interpretation.The derived results are of course equally relevant for CT with other radiation sources of low brilliance, notably neutron tomography.

FIG. 3 .
FIG.3.Exponential and linear convolution.The example structure h is convolved (by * logexp and * ) with a normalized uniform density of width 2s x .The lower half shows the relative deviation a = (j * logexp ah − j * ah)/a for different absorption strengths a.The numbers a = 1, 3, and 5 correspond to approximately 37%, 5%, and 0.7% transmission, respectively.

g
blur := M −1 L s y −s y S ys j (r s )p(r s ,r d )dx s dy s dx d = s y −s y j (y s ) M −1 L S ys j y s (x s )p(r s ,r d )dx s dx d dy s = s y −s y j (y s )g w (n(y s ),s(y s ))dy s .(25) Note that y s → n(y s ) and y s → s(y s ) are smooth and for small cone angles can be approximated by n(y s ) ≈ n(0) and s(y s ) ≈ s(0) + (M − 1)M −1 y s .Defining jeff (y s ) := M(M − 1) −1 j (−M(M − 1) −1 y s ) analogously to j eff , Eq. (25) can be rewritten as s )g w (n(0),s(0) − y s )dy s (26) = [ jeff * g w (n(0),•)](s(0)), 3 denote the voxels of the reconstruction volume.The fully discrete version of Eq. (29) then becomes

FIG. 6 .
FIG. 6. Experimental setup and results.(a) Photograph showing the setup consisting of the sealed tube source, kappa goniometer, and pixel detector.The background has been masked to highlight the components.(b) "Raw" projection prior to the reconstruction.The bar at the top indicates the effective width of the source 2s x .(c) Numerical projection of the reconstructed volume.(d) Slice through the reconstructed volume.(e) Magnified region [outlined in (d)].
(b)], the 3DRT reconstruction yields a very sharp reprojection of the reconstruction [Fig.6(c)] and shows surprising interior detail [Figs.6(d) and 6(e)].Effects of the nonlinear error in the form of some cloudy background around the edges of strongly absorbing bones can be identified; they, however, do not significantly hamper the overall reconstruction quality.

TABLE I .
Geometric parameters of the tomography setup: source to sample distance z 01 , sample to detector distance z 12 , geometric magnification M, cone angle σ , effective source extents (in the object plane) s x ,s y , and effective pixel size (in the object plane) p eff .