A GENERAL SCHEME FOR CONSTRUCTING INVERSION ALGORITHMS FOR CONE BEAM CT

Given a rather general weight function n0, we derive a new cone beam transform inversion formula. The derivation is explicitly based on Grangeat’s formula (1990) and the classical 3D Radon transform inversion. The new formula is theoretically exact and is represented by a 2D integral. We show that if the source trajectory C is complete in the sense of Tuy (1983) (and satisfies two other very mild assumptions), then substituting the simplest weight n0 ≡ 1 gives a convolution-based FBP algorithm. However, this easy choice is not always optimal from the point of view of practical applications. The weight n0 ≡ 1 works well for closed trajectories, but the resulting algorithm does not solve the long object problem if C is not closed. In the latter case one has to use the flexibility in choosing n0 and find the weight that gives an inversion formula with the desired properties. We show how this can be done for spiral CT. It turns out that the two inversion algorithms for spiral CT proposed earlier by the author are particular cases of the new formula. For general trajectories the choice of weight should be done on a case-by-case basis.


Introduction.
The inversion of cone beam transform in R 3 is interesting both from a theoretical point of view (as a problem in integral geometry) and a practical point of view (as a problem in computed tomography).Many important results are known in this area.A very incomplete list includes microlocal analysis of the inversion problem [8,16], relationships between the cone beam and Radon transforms [6,7,15,17,18], completeness condition [20], and inversion formulas [5,6,20].The introduction of spiral scanning and two-dimensional detector arrays brought about a great number of works devoted to various approximate and exact reconstruction algorithms for spiral CT (see, e.g., review paper [19] and the references therein).The major problem in spiral CT is that until recently, there did not exist an inversion formula that admitted shift-invariant filtering and back projection (FBP) implementation.Following the approach of [3,13], several FBP-type inversion algorithms have been proposed in the literature.However, they were either approximate or used shift-variant filtering.The computational efficiency is very important in practice because contemporary scanners generate huge amount of data and inefficient algorithms are not capable of reconstructing the image in a reasonable time.
In [10,11,12], two new formulas for inverting the cone beam transform in the case of spiral source trajectory have been proposed.An important feature of these formulas is that they are theoretically exact, have an FBP structure, and the filtering step is shift-invariant.Even though the two inversion formulas are proven to be exact, the proofs contained in [10,12] are not very illuminating.First, they do not show where the formulas come from.Second, even though it is intuitively clear that the results of [10,11,12] should somehow be related to the Radon transform in R 3 , the presented proofs do not allow to see that.Third, it is not clear how to generalize these results to other source trajectories.
In this paper we derive a more general cone beam inversion formula.The derivation is explicitly based on (a) Grangeat's relation [7] between the cone beam and Radon transforms in R 3 and (b) the classical 3D Radon transform inversion.The distinctive features of our approach can be summarized as follows: (1) the choice of a more general weight n 0 (s,x,α); (2) a novel way of dealing with discontinuities in the weight; (3) a novel way of eliminating the intermediate function (i.e., the first derivative of the Radon transform).Given a Radon plane Π, one usually assigns weights to the points of the intersection of Π with the source trajectory C independently of which point x ∈ Π is currently reconstructed (see, e.g., [14]).Here, to the contrary we assume that given x ∈ Π, the distribution of weights among the points of intersection depends on x.
Any approach to use the 3D Radon transform theory to obtain an FBP type cone beam transform inversion formula deals with the problem of discontinuities of the weight in one way or another.For example, the idea of [3] is to replace a discontinuous weight function by a smoother one, which still guarantees theoretical exactness.Here, we use a different approach.At the beginning, we introduce a smooth cutoff function η that eliminates all the discontinuities and singularities (resulting, e.g., from the points of tangency).Then, at the very end, we consider the limit as η → 1.It turns out that these discontinuities are precisely what shows up in the inversion formula.
The method of eliminating the intermediate function is based on a derivation analogous to the one in [9].The key observation is that the derivative of the cone beam transform, which appears in Grangeat's formula, can be removed from that formula and applied to a factor φ appearing in front of it.It turns out that under certain assumptions about the weight, φ is piecewise constant.So, in the limit as η → 1, one of the integrals reduces to a finite sum over the discontinuities of φ.This explains why the new inversion formula is represented by a 2D integral, whereas we started with a 3D one (one integration in Grangeat's relation and a two-dimensional integral in the 3D Radon transform inversion formula).
Note that the inversion formula is derived under some rather general assumptions about the weight n 0 .By construction, the formula is theoretically exact and is represented by a 2D integral.At the outset, it neither has the FBP structure nor the filtering step is shift-invariant.We show that if the source trajectory C is complete in the sense of Tuy [20] (and satisfies two other very mild assumptions), then substituting the simplest weight n 0 ≡ 1 gives a convolution-based FBP algorithm.However, this easy choice is not always optimal from the point of view of practical applications.The weight n 0 ≡ 1 works well for closed trajectories, but the resulting algorithm does not solve the long object problem if C is not closed.In the latter case one has to use the flexibility in choosing n 0 and find the weight that gives an inversion formula with the desired properties.We think that it is probably impossible to give a recipe for finding the weight that gives an optimal convolution-based FBP algorithm for an arbitrary complete trajectory C. Instead, this should be done on the case-by-case basis.We show how this can be done for spiral CT.It turns out that the results of [10,12] are simply two particular cases of the new formula.
In Section 2, we derive the new inversion formula.In Section 3, we show that for a general C the weight n 0 ≡ 1 yields a convolution-based FBP algorithm.As an example we consider the classical two-orthogonal-circles trajectory.In Section 4, we consider two particular weights for spiral CT that yield convolution-based FBP algorithms.The latter coincide with the algorithms of [10,12].
Throughout the paper, it is tacitly assumed that all the projections required for the algorithm are available.The flexibility in choosing n 0 can also be used for reducing the amount of projections needed for the algorithm.For example, both weights for spiral CT, which are discussed here, allow the use of axially truncated projections.Additionally, one weight requires a detector array smaller than the other.
Introduce the sets Crit(s, x) An important ingredient in the construction of the inversion formula is weight function Remark 2.4.The function n 0 can be understood as follows: x and α determine the plane Π(x, α), and the weight n 0 assigned to y(s) ∈ Π(x, α) ∩ C depends on the location of x (see Figure 2.1).In view of this interpretation, we assume that n 0 (s,x,α) = n 0 (s, x, −α). Define In (2.4) and throughout the paper, j denotes the sum over all s j such that The main assumptions about n 0 are the following properties.
Property 2.6.There exist finitely many The inversion formula, which is to be derived here, holds pointwise.Therefore, if f needs to be reconstructed for all x belonging to a set U , then Properties 2.1, 2.2, 2.3, 2.5, and 2.6, are supposed to hold pointwise, and not uniformly with respect to x ∈ U.
Alternatively, (2.16) can be derived as follows.Denote, temporarily, (2.17) Then, the integral with respect to α in (2.11) can be written as follows: where dθ is the Lebesgue measure on the great circle β ⊥ (s, x) ⊂ S 2 .Clearly, (2.18) is equivalent to (2.16).In the step from line 1 to line 2 in (2.18), we have tacitly extended A(s, x, α) to a small neighborhood of {(s, α) : s ∈ I, α ∈ β ⊥ (s, x)}.Since η(α) erases all the singularities and n(s, x, α) is locally a constant in a neighborhood of all regular points, this can be easily done.The step from line 2 to line 3, which involves changing the order of integration, can be justified by noticing that A(s, x, α) ≡ 0 in a neighborhood of all α such that α • ẏ(s) = 0 for some s ∈ I and using the Lebesgue dominated convergence theorem.
Because of Properties 2.2 and 2.6, the number of terms in the summation in (2.28) is finite.Due to the factor sgn(α • ẏ(s)) in (2.27), this number is at most one plus the number of discontinuities of n(s, x, α).Also, c m (s, x) = 0 if s ∈ I \ I reg (x) (i.e., Crit(s, x) = β ⊥ (s, x)).This follows from construction, because η(α) ≡ 0 on β ⊥ (s, x) for such a pair (s, x).
3. The case n 0 ≡ 1.In this section, we suppose that n 0 (s,x,α) ≡ 1 and show that this choice of n 0 leads to a convolution-based FBP algorithm.First, verify Properties 2.5 and 2.6.Since n Σ (x, α) is just the number of points in Π(x, α)∩ C, assuming that Property 2.1 holds for all x ∈ supp f , we immediately get n Σ (x, α) ≥ 1 on suppf × S 2 .From (2.4), n(s, x, α) = 1/n Σ (x, α).Taking s ∈ I reg (x), and restricting α to β ⊥ (s, x), we see that this ratio can be discontinuous only if Π(x, α) is tangent to C or contains an endpoint of C. As in the previous section, such a plane is called critical.Property 2.6 follows from Property 2.2.
In addition, if α ∈ β ⊥ (s, x), the ratio n(s, x, α) = 1/n Σ (x, α) depends on x only via β(s, x).The same, clearly, applies to φ(s, x, α) defined in (2.27).Consequently, c m (s, x) and α ⊥ (s,x,θ m ) depend on x only via β(s, x).Therefore, we can replace x by β(s, x) in the arguments of c m and α ⊥ and rewrite (2.28) in the form At this point, however, it is not clear that the filtering step defined by (3.2) is shift-invariant.Now, fix y(s) ∈ C, and let Π 0 y(s) be a critical plane for some x 0 ∈ supp f .It is obvious that all other x ∈ Π 0 will share this plane as critical.Furthermore, the set of all critical planes Π 0 y(s) can be described as a finite union of oneparametric families of planes.Each endpoint of C generates a family, and each smooth segment of C generates a family of tangent planes.
Again, fix a critical plane Π 0 y(s) and pick any x ∈ Π 0 .Then, for the appropriate m, the unit vector α ⊥ (s, β(s, x), θ m ) is parallel to Π 0 (cf.Remark 2.7).By construction, the vectors cos γβ(s, x) + sin γα ⊥ (s, β(s, x), θ m ), 0 ≤ γ < 2π , belong to the same plane Π 0 .Let ω be the polar angle in that plane.Since we can write (with abuse of notation) Therefore, The integral in (3.4) is a convolution.Therefore, one application of FFT to (3.4) gives values of Ψ m (s, β(s, x)) for all x ∈ Π 0 at once.Equations (3.1) and (3.4) imply that the resulting algorithm is of the FBP type.The first step is to perform convolution filtering of the derivative of cone beam projections according to (3.4) for all required m and Π 0 .The second step is to backproject the filtered data according to (3.1).
As an example, consider how this theory applies to the classical trajectory consisting of two orthogonal circles.In this case, where supp f is supposed to be inside a ball centered at the origin and with sufficiently small radius r < R. The (virtual) detector plane moves simultaneously with the source and, at each instant, its equation is {z ∈ R 3 : z • y(s) = 0}.Since the cases y(s) ∈ C 1 and y(s) ∈ C 2 are completely analogous, we consider only the case y(s) ∈ C 1 .From (2.27) with n = 1/n Σ , φ is discontinuous if Π 0 y(s) is parallel to ẏ(s) or is tangent to C 2 .In the first case, the magnitude of the jump is 1/2.Indeed, such a plane has four points of intersection with C, so φ = 1/4 on the side of the jump, where α • ẏ(s) > 0 and φ = −1/4 on the other side.The corresponding direction of integration is parallel to ẏ(s) (see the dot-dashed lines in Figure 3.1).In the second case, the number of intersections changes from two to four, so the magnitude of the jump is 1/4.The corresponding set of integration directions is illustrated in compute the contribution from the endpoints to the image.If the trajectory C is short (e.g., as in C-arm scanning), this should not cause any problems.However, if C is long (e.g., as is frequently the case in spiral scanning), using the endpoints for image reconstruction at any given x inside the ROI leads to undesirable consequences: long-object problem is not solved, requirements on the detector array are excessive, and others.This argument shows that for long-source trajectories image reconstruction at x should be performed using a section C(x) ⊂ C.However, this brings additional difficulties.If the entire ROI is split into a small number of sub-ROIs and C(x) is the same inside each sub-ROI, artifacts are likely to appear at the boundaries of the sub-ROIs.On the other hand, if C(x) changes smoothly with x, then (3.1) and (3.2) do not apply, because the dependence of φ(s, x, α) on x is no longer via only β(s, x).
Our discussion shows that for long trajectories, the weight n 0 ≡ 1 is not optimal.The way out is to find a more complicated (piecewise constant) weight that would cancel out the contributions from the endpoints.In the next section, we demonstrate these weights for spiral CT.

Construction of weights for spiral CT. Consider a spiral path of the Xray source
As was shown in [2,4], any point strictly inside the spiral belongs to one and only one PI segment.Recall that a PI segment is a segment of line endpoints of which are located on the spiral and separated by less than one pitch in the axial direction.Let s = s b (x) and s = s t (x) denote values of the parameter corresponding to the endpoints of the PI segment containing x.We call I PI (x) := [s b (x), s t (x)] the PI parametric interval.The part of the spiral corresponding to I PI (x) is denoted by C PI (x).As C(x), which is used for image reconstruction at x, we take the spiral segment C PI (x).It is clear that any plane through x intersects C PI (x) at least at one point.Also, inside the PI parametric interval there exists š = š(x) such that the plane through y(š) and parallel to ẏ(š) and ÿ(š) contains x.It is assumed throughout this section that f ∈ C ∞ 0 (U ), where , is an open cylinder located strictly inside the spiral.
Since C PI (x) depends on x, (2.28) and (2.29) apply if C PI (x) has Properties 2.1, 2.2, and 2.3 for every x ∈ U .As was already mentioned, C PI (x) is complete for any x ∈ U .The other two properties are quite obvious as well.[12].The weight n 0 of this subsection depends on two auxiliary vector functions e k (s, β), k = 1, 2. The first one is easy;

Derivation of inversion formula (2.13) of
2) e 1 (s, β) ∈ β ⊥ is a unit vector in the plane through y(s) and spanned by β and ẏ(s).
The second vector function is more complicated.Given y(s), s ∈ (s b (x), s t (x)) \ {š(x)}, find s tan ∈ I PI (x) and s tan ≠ s such that the plane through x, y(s), and y(s tan ) is tangent to C PI (x) at y(s tan ).Existence and uniqueness of such s tan = s tan (s, x) are shown in [12].Once s tan has been found, e 2 (s, x) ∈ β ⊥ (s, x) is a unit vector in the plane through x, y(s), and y(s tan ).Vector function e 2 (s, x) can be extended with respect to s to all of I PI (x) as a continuous function.The direction of e 2 (s, x) is chosen so that e 1 (s, x) = e 2 (s, x), when s = š(x) (see [12]).
Fix s ∈ I PI (x).Pick any x ∈ U such that s ∈ I PI (x), and find e 2 = e 2 (s, x).It can be shown that the same vector e 2 works for all x ∈ U with β(s, x) = β(s, x ).In particular, s ∈ I PI (x ) for all such x .Thus, we have actually determined e 2 (s, β), where β = β(s, x), s ∈ I PI (x).Define It is proven in [12]  The integral with respect to γ in (2.28) is odd when θ → θ +π .Similarly, the values of the jump of φ at two points θ m 1 and θ m 2 separated by π differ by a factor −1. So, by inserting an extra factor 2, this integral can be confined to an interval of length π .This implies that we can take α ⊥ (s,β,θ m ) = e m (s, β), m = 1, 2, and (2.28) transforms into the inversion formula of [12] (note that all jumps of φ have amplitude 1): where

Derivation of inversion formula (2.13
) of [10].The weight n 0 of this subsection depends on only one auxiliary vector function e(s, β).Choose any ψ ∈ C ∞ ([0, 2π]) with the properties Suppose that s 0 , s 1 , and s 2 are related by Since ψ(0) = 0, s 1 = s 1 (s 0 ,s 2 ) is a continuous function of s 0 and s 2 .Conditions (4.7) and (4.8) imply that s 1 ≠ s 2 unless s 0 = s 1 = s 2 .In order to avoid unnecessary complications, we assume also that ψ (0) = 0.5; ψ (2k+1) (0) = 0, k≥ 1. (4.9) If (4.9) holds, then s 1 = s 1 (s 0 ,s 2 ) is a C ∞ -function of s 0 and s 2 .Conditions (4.7) and (4.9) are very easy to satisfy.We can take, for example, ψ(t) = t/2, and this leads to It is shown in [10] that u(s 0 ,s 2 ) is a C ∞ -vector function of its arguments.Fix x ∈ U and s 0 ∈ I PI (x).Find s 2 ∈ I PI (x) such that the plane through y(s 0 ), y(s 2 ), and y(s 1 (s 0 ,s 2 )) contains x.More precisely, we have to solve for s 2 the following equation: x − y s 0 • u s 0 ,s 2 = 0, s 2 ∈ I PI (x).(4.12) It is shown in [10] that such s 2 exists, is unique and depends smoothly on s 0 .Again, it can be seen that if e(s 0 ,x), s 0 ∈ I PI (x) is found for x ∈ U , then the same vector works for all x ∈ U with β(s 0 ,x) = β(s 0 ,x ) (see [10]).Therefore, this construction defines s 2 := s 2 (s 0 ,β) and, consequently, u(s 0 ,β) := u(s 0 ,s 2 (s 0 ,β)).Finally, we set e(s, β) It is proven in [10]  Arguing in the same way as before, we immediately obtain the inversion formula of [10] f  With the choice of the weight n 0 given by (4.3) or (4.13), φ(s, β(s, x), θ) is continuous in θ as the planes Π(y(s), α), α = α(s, θ) ∈ β ⊥ (s, x) cross the boundary of C PI (x).Therefore, no contribution from the endpoints of C PI (x) to the image at x needs to be computed.Alternatively, we can say that the contributions from the endpoints to the image cancel each other.
The two examples presented in this section show that it is important to have the extra flexibility provided by the dependence of n 0 on x.Indeed, as is seen from (4.3) and (4.13), n 0 (s,x,α) explicitly depends on x.
Since the weight n 0 in the general inversion formulas (2.28) and (2.29) is arbitrary (subject to Properties 2.5 and 2.6), we can use this flexibility to improve the detector usage.For example, (4.15) provides much better detector usage than (4.5) and (4.6) (see [10]).

Call for Papers
Thinking about nonlinearity in engineering areas, up to the 70s, was focused on intentionally built nonlinear parts in order to improve the operational characteristics of a device or system.Keying, saturation, hysteretic phenomena, and dead zones were added to existing devices increasing their behavior diversity and precision.In this context, an intrinsic nonlinearity was treated just as a linear approximation, around equilibrium points.
Inspired on the rediscovering of the richness of nonlinear and chaotic phenomena, engineers started using analytical tools from "Qualitative Theory of Differential Equations," allowing more precise analysis and synthesis, in order to produce new vital products and services.Bifurcation theory, dynamical systems and chaos started to be part of the mandatory set of tools for design engineers.
This proposed special edition of the Mathematical Problems in Engineering aims to provide a picture of the importance of the bifurcation theory, relating it with nonlinear and chaotic dynamics for natural and engineered systems.Ideas of how this dynamics can be captured through precisely tailored real and numerical experiments and understanding by the combination of specific tools that associate dynamical system theory and geometric tools in a very clever, sophisticated, and at the same time simple and unique analytical environment are the subject of this issue, allowing new methods to design high-precision devices and equipment.
Authors should follow the Mathematical Problems in Engineering manuscript format described at http://www .hindawi.com/journals/mpe/.Prospective authors should submit an electronic copy of their complete manuscript through the journal Manuscript Tracking System at http:// mts.hindawi.com/according to the following timetable:

First
Round of ReviewsMarch 1, 2009 (restricted to s ∈ I PI (x)) is discontinuous, when α(s, θ) is perpendicular to either e 1 (s, β(s, x)) or e 2 (s, β(s, x)).Equation (4.4) and the argument preceding (4.3) imply that if s ∈ I PI (x), then c m (s, x) and α ⊥ (s,x,θ m ) depend on x only via β(s, x).In fact, whether the inclusion s ∈ I PI (x) itself holds or not is determined only by β(s, x).So, (3.1) and (3.2) do apply when I is replaced by I PI (x).