A reconstruction approach for imaging in 3D cone beam vector field tomography.

3D cone beam vector field tomography (VFT) aims for reconstructing and visualizing the velocity field of a moving fluid by measuring line integrals of projections of the vector field. The data are obtained by ultrasound measurements along a scanning curve which surrounds the object. From a mathematical point of view, we have to deal with the inversion of the vectorial cone beam transform. Since the vectorial cone beam transform of any gradient vector field with compact support is identically equal to zero, we can only hope to reconstruct the solenoidal part of an arbitrary vector field. In this paper we will at first summarize important properties of the cone beam transform for three-dimensional solenoidal vector fields and then propose a solution approach based on the method of approximate inverse. In this context, we intensively make use of results from scalar 3D computerized tomography. The findings presented in the paper will continuously be illustrated by pictures from first numerical experiments done with exact, simulated data.

1. INTRODUCTION

Vector field tomography (VFT) deals with the problem of reconstructing a vector field, for example, a velocity field of an incompressible, moving fluid, from line integrals of projections of the field. VFT has various applications in photoelasticity, oceanography, nondestructive testing, and medical imaging, where we may think of tumor detection by reconstructing and visualizing blood flow which is known to be more irregular and more intense around tumors than in normal tissue, see [1]. The integral data can be measured using ultrasound signals when we assume that the Doppler shift of the frequency is approximately proportional to the velocity of the particle in the fluid which causes the shift. This is a reasonable assumption if the particle velocity is significantly smaller than the speed of sound within the medium under consideration.

Although this seems to be quite simple at first sight it will become clear after the definition of the cone beam transform for vector fields that only the projection of the vector field onto the line of integration can be measured. This enormous loss of information is the reason why we can only hope to recover the solenoidal part of the vector field from our measurements, a fact which has for example been shown in [2].

Ultrasound devices may be utilized as a supportive method in preliminary examinations for tumor detection by reconstructing and visualizing blood flow which has already been suggested in 1977 by Wells et al. [3]. This may help to reduce the radiation dose of a patient tremendously. Thinking of mammography, it is known that the pressure which is put on the breast during the examination may crush small existing tumors and allow them to spread more easily. This danger could also be avoided or at least reduced by using ultrasound for preventive medical examinations. It should be noted that these deliberations are subject to the condition that the algorithms and the medical equipment based on ultrasound work as reliable and fast as current X-ray techniques.

A possible measurement setup where the scanning curve Γ is a circle in the plane {x 3 = 0} is depicted in Figure 1. This defines exactly the geometry used in our first numerical experiments which are presented in the last section of this work. A lot of theoretical and numerical results have been achieved over the last few years for the parallel geometry. Juhlin [4] suggested a measurement setup which is suited to fully reconstruct solenoidal fields in two dimensions. Mathematical properties of this model can be found in Sparr et al. [5]. The singular value decomposition for the 2D fan-beam Radon transform of tensor fields has been presented in an article by Kazantsev and Bukhgeim [6]. Desbat and Wernsdörfer [7] developed an iterative method. For 3D Doppler tomography, Schuster established an inversion scheme of filtered backprojection type [8, 9] relying on the method of approximate inverse. Together with Rieder [10], he obtained convergence with rates and stability with respect to noisy data for this method.

Figure 1

Measurement setup using a circle in the plane {x 3 = 0} as scanning curve Γ.

As in scalar 3D computerized tomography, the cone beam transform is of special interest from a practical point of view. It is defined for a tensor field of rank m by where α ∈ Γ is a source point on the scanning curve which surrounds the object Ω, ω ∈ S is the unit vector of direction of the line and f is a tensor field of rank m with compact support in the open domain Ω. Tensor fields of rank m = 0 are scalar functions f(x), tensor fields of rank m = 1 are vector fields f(x) in ℝ. In (1) we use Einstein's summation rule, that means we sum up over equal indices i , where 1 ≤ i ≤ n.

A lot of theory has been developed for the common case of tensor fields of arbitrary rank. To facilitate notation and to direct the readers' attention to the cases important for practical applications we confine the following remarks to the scalar cone beam transform as well as the cone beam transform for vector fields in n dimensions. In the following, scalar fields will be denoted by f whereas vector fields will be written bold-faced, such as f.

Setting m = 0 in (1) we obtain the well-known (scalar) cone beam transform of a scalar field f : ℝ ⊃ Ω → ℝ. For m = 1 formula (1) is the cone beam transform for vector fields f : ℝ ⊃ Ω → ℝ, which in 3D is also often referred to as the Doppler transform. It reads as Hence the mathematical problem of 3D cone beam VFT consists of inverting D 1 f = g for given measurements g ∈ ℝ, that is reconstructing a three-dimensional vector field f from one-dimensional, that is, scalar, data g. In contrast to D 0 an inversion formula for D 1 is not known by now and an inversion scheme for the cone beam transform for vector fields has not been established so far.

From (3) it can be seen that only the integral of the projections of the vector field along the ray of integration can be obtained from the measurements, which is emphasized in Figure 2. Considering an ultrasound wave starting from the source point α on the scanning curve Γ only the projection of the green sample vector onto the line of integration can be measured. The projection is illustrated by the red arrow. This also means that any vector orthogonal to the ultrasound wave does not contribute to the integral at all, a property which is depicted by the unit circle and the corresponding sample vectors orthogonal to ω. As a consequence, a full reconstruction of an arbitrary vector field is impossible with the underlying measurement geometry. As already mentioned in the abstract of our paper, we can only hope to reconstruct the solenoidal part of an arbitrary vector field since the vectorial cone beam transform of any gradient vector field with compact support is identically equal to zero. This can easily be seen from the following short computation for a gradient field f(x) = ∇ϕ(x) = (∂ ϕ(x),…, ∂ ϕ(x))⊤, ϕ : ℝ ⊃ Ω → ℝ, with compact support

Figure 2

Visualization of the measured projection data for D 1 f.

The method of approximate inverse introduced by Louis and Maass [11] delivers a mathematical framework for coping with inverse problems in an efficient way. The method computes a smoothed version of the solution f with the help of so-called mollifiers. These are smooth approximations to delta functions. Using a duality argument the method then consists of evaluations of inner products of the given data g with reconstruction kernels. One feature of the method is, that invariances of the underlying operator can be used to speed up the computation time tremendously. This method was successfully applied to the reconstruction problem in 3D computerized tomography, that is, (2), see [12]. We use these results to extend the method of approximate inverse to (3).

We summarize the contents of the paper. First we state essential mathematical properties of (1) for the special cases of m = 0 and m = 1. The interested reader will find generalizations of the presented theorems to symmetric, covariant tensor fields of any rank and in any dimension in [13], especially the extension of Grangeat's formula which has been proven for (1) in case n = 3, m = 0 (see (2)) in [14]. Then we outline how the method of approximate inverse can be used for (2) to solve D 0 f = g and present its application to (3) to solve D 1 f = g. Furthermore, an approach is presented, how reconstruction kernels for D 1 from (3) can be calculated with the help of known reconstruction kernels for the scalar cone beam transform D 0. Some pictures from numerical experiments show that this approach is very promising.

2. MATHEMATICAL PROPERTIES OF D 0 AND D 1

Let denote the space of square integrable functions from X ⊂ ℝ to ℝ where the L 2-inner product of two functions is given as Fundamental properties of (2) and (3) are summarized in the following theorem, which is the result of straightforward calculations.

Theorem 1

Let Ω := {x ∈ ℝ : |x| < 1} with ∂Ω = S the reconstruction region, that is, the region in which the object Ω is contained. The mappings are linear and bounded if The adjoints ( D 0 ∗ : L 2(Γ × S ) → L 2(Ω) and D 1 ∗ : L 2(Γ × S ) → L 2(Ω, ℝ) are given by

For n = 3, we obtain the well-known cone beam transform D 0 with the corresponding backprojection operator D 0 ∗ of scalar fields which is thoroughly investigated in 3D computerized tomography as well as the mathematical model of 3D cone beam vector tomography and the backprojection reads as D 1 ∗ g(x) represents an integration over all lines intersecting x. The result can be seen as an average over all directions connecting a source point α and x.

One of the crucial tools when computing reconstruction kernels in scalar cone beam tomography is the formula of Grangeat [14]. We proved a generalization of that formula which is valid for any symmetric, covariant tensor field of rank m in n dimensions in [13] but our presentation will be restricted to D 0 and D 1.

Theorem 2 (Schuster [13] based on Hamaker et al. [15]).

Assume n ≥ 2, f ∈ 𝒞 0 ((Ω), and f ∈ 𝒞 0 ((Ω, ℝ). Then, where α ∈ Γ, ω ∈ S , dS denotes the surface measure on S , R is the and is the projection of f onto (x − α)/ |x − α|.

In the following, we want to limit our presentation to n = 3. Then formula (12) reads as

In the scalar case a solver for D 0 can be constructed with the help of (11) for n = 3. This is done by inverting the Radon transform R which is possible if the condition of Tuy-Kirillov is satisfied. It tells that we have full knowledge of R g(ω, s) for all ω, s and any scalar function g : Ω3 → ℝ, if any plane intersecting the object Ω ⊆ Ω3 does also have at least one intersection point with the scanning curve Γ and this intersection must be nontransversally. This works fine for D 0 since then f(x) is independent of α but unfortunately that does not help in case of D 1 (and analogous transforms for tensor fields as well, see [13]), since there the object function f = 〈f(x), (x − α)/|x − α|〉ℝ of R depends on and hence changes with α, see (15). Thus we seek an alternative way of solving D 1 f = g for vector fields f.

3. APPROXIMATION OF RECONSTRUCTION KERNELS IN VECTOR FIELD TOMOGRAPHY

As already said, an inversion formula for D 1 is not known by now and an inversion scheme for the cone beam transform for vector fields has not been established so far. So, the aim of this paper is to deduce a completely new method for three-dimensional cone beam VFT. A comparison with existing algorithms is difficult since there are no inversion methods for the vectorial cone beam transform to the authors' best knowledge. Nevertheless, some famous algorithms will certainly come to the reader's mind when thinking of tomography. The well-known FDK algorithm (Feldkamp et al. see [16]) is detailed on [17, page 128]. From this description it immediately becomes clear that the algorithm does not work for D 1. The integrand of the cone beam transform for vector fields strongly depends on the direction ω, a fact which is explicitly disregarded by the FDK algorithm. The methods of Norton (see [18, 19]) and Prince [20] are specifically suited to solve 2D, respectively, 3D problems for vector fields in parallel geometry. They both use transforms different from D 1. The generalization of Norton's approach to 3D vector tomography of Doppler-transformed fields in parallel geometry was a challenging problem for Lade et al. in [21]. Regrettably, neither approach can be adapted to VFT using the cone beam geometry. Finally, no Fourier slice theorem for VFT is known, not even for standard 3D cone beam tomography, so Fourier methods are not an alternative.

The method of approximate inverse, which was established by Louis and Maass [11] in 1990, leads to an algorithm of filtered backprojection type if invariances and some appropriate approximations are used. This has been shown for example in [22, 23] or [12]. Fundamental properties of it have also been published in [24, 25]. Its theory was enhanced over the last decade and the method was successfully applied to several reconstruction problems in medical imaging and nondestructive testing, such as computerized tomography, inverse scattering, thermoacoustic computerized tomography, diffractometry, and Doppler tomography. In [12] the method was applied to 3D cone beam tomography, that is, to D 0. We describe this approach and then formulate an extension of it to D 1.

Let f ∈ L 2(Ω3) be a scalar function. The approximate inverse computes a smoothed version f of f by convolving f with a mollifier e ∈ 𝒞 ∞(ℝ3). A mollifier e is a smooth function with small essential support having the property that Here, ∗ denotes the convolution Such a function is given by the Gaussian kernel Provided that we can solve the equation then we can reconstruct f from the measured cone beam data D 0 f by where v (x) = v (x; α, ω) ∈ L 2(Γ × S 2) for x ∈ Ω is called a reconstruction kernel. Hence the method of approximate inverse consists of evaluating inner products of the given data D 0 f with reconstruction kernels v (x). This can be done in an efficient way using the translation invariance of e and some appropriate approximations which are outlined in detail in [12]. These imply that we have to solve (19) only once, namely for x = 0, and apply the invariances to get the remaining reconstruction kernels. By doing so the computation time is shortened significantly. The computation of reconstruction kernels for circular 3D cone beam tomography, that is, for exactly the same scanning geometry that we use for the reconstruction of vector fields, has been detailed in [26].

In comparison to other regularization methods such as Tikhonov regularization which would result in enormous computational costs because of the very large, full matrices, the approximate inverse is much more efficient, an advantage that is especially crucial in tomographic applications because of the large number of evaluations.

To apply the method to D 1 and hence to VFT, we construct mollifier fields E ∈ L 2(Ω3, ℝ3) defining where e 1 = (1, 0, 0)⊤, e 2 = (0, 1, 0)⊤ and e 3 = (0, 0, 1)⊤. Using again the Gaussian (18) as mollifier e we obtain for f ∈ L 2(Ω3, ℝ3). Unfortunately, by now the exact reconstruction kernels V (x), that is, the solutions of D 1 ∗[V (x)] = E (x − ·) are still unknown. But the special structure of the mollifier fields E allow for a computation of reconstruction kernels for with the help of kernels for D 0.

Theorem 3

Let v be the reconstruction kernel associated to e with respect to D 0, that is, Defining V (x; α, ω) := v (x; α, ω)·e ∈ L 2(Γ × S 2, ℝ3) yields that means V is a reconstruction kernel associated to E with respect to P. The adjoint P ∗ of P is given as for g ∈ L 2(Γ × S 2, ℝ3).

Proof

The adjoint P ∗ is computed as where we applied Fubini's theorem as well as the substitution x = α + t ω. A short calculation further shows that

The data Pf are not known and cannot be computed from D 1 f. But, observing that and since Pf(α, ω) ∈ ℝ3 we obtain where ω 1 ⊥, ω 2 ⊥ ∈ S 2 are such that {ω, ω 1 ⊥, ω 2 ⊥} is an orthonormal basis of ℝ3 and λ 1, λ 2 are appropriate coefficients. Thus approximating we neglect the parts orthogonal to ω and can apply the method of approximate inverse using the reconstruction kernels V for P. This procedure results in Algorithm 1.

Algorithm 1

For cone beam VFT.

It is worth to mention that the mathematical model has not been changed. Results for the scalar cone beam transform D 0 f(α, ω) are transferred to its 3D equivalent Pf(α, ω) which is an approximation to D 1 f(α, ω) by using Pf(α, ω) ≈ D 1 f(α, ω)ω.

Figures 3 and 4 display first results of the above algorithm when applied to exact, simulated data for the solenoidal vector fields f(x) = (−x 2, x 1, 0)⊤ and f(x) = (1 − x 2 2 − x 3 2, 0, 0)⊤. As already said, we made use of the measurement setup as shown in Figure 1. The scanning curve was Γ = r S 2 ∩ {x 3 = 0}, r = 3, that is a circle of radius r = 3 in the plane {x 3 = 0}. The divergence-free vector fields are assumed to be defined in the three-dimensional unit ball, that is, according to our definitions we have Ω = Ω3. The reconstructions depicted in Figures 3 and 4 were done in the plane {x 3 = 0}. The mollifier e defining the fields E was chosen as the Gaussian given in (18). The regularization parameter was γ = 0.00692 and γ = 0.007, respectively. The corresponding reconstruction kernel v (0; ·, ·) is depicted in Figure 5. The reconstruction kernel can be seen as a lowpass filter. Then, the regularization parameter γ determines the width of the filter and thus can be interpreted as the cutoff frequency. Large values of γ correspond to a large smoothing effect in the reconstructed vector field. Unfortunately we cannot estimate the range of an optimal γ since it depends on the noise level of the measured data as well as on the exact solution f itself. Nevertheless, in our experiments a value of γ ≈ 0.007 always led to good results.

Figure 3

(a) Original vector field f(x) = (−x 2, x 1, 0)⊤. (b) Reconstruction with the described algorithm for γ = 0.00692 using exact, simulated data D 1 f.

Figure 4

(a) Original vector field f(x) = (1 − x 2 2 − x 3 2, 0, 0)⊤. (b) Reconstruction with the described algorithm for γ = 0.007 using exact, simulated data D 1 f.

Figure 5

Reconstruction kernel v (0; α, ω) for γ = 0.007 associated to the scalar cone beam transform D 0 according to [12].

Figure 4 emphasizes that the main part of the reconstruction error is located at the boundary of the domain. Although our scanning curve is only one circle our algorithm nevertheless allows us to reconstruct any arbitrary plane in the x 3-direction. Figures 6 and 7 show the reconstructions for the planes {x 3 = 0.5}, {x 3 = 0.75}, {x 3 = 0.9} and {x 3 = 0.95} of the two afore-mentioned vector fields. Figure 7 also illustrates that the intensity of the vector field f(x) = (1 − x 2 2 − x 3 2, 0, 0)⊤ decreases as should be expected since x 3 2 is subtracted in the first component and that even the directional error at the boundary is reduced with increasing x 3.

Figure 6

Reconstruction of the vector field f(x) = (−x 2, x 1, 0)⊤ with the described algorithm for γ = 0.00692 using exact, simulated data D 1 f. (a),(b) The planes {x 3 = 0.5} and {x 3 = 0.75}. (c),(d) The planes {x 3 = 0.9} and {x 3 = 0.95}.

Figure 7

Reconstruction of the vector field f(x) = (1 − x 2 2 − x 3 2, 0, 0)⊤ with the described algorithm for γ = 0.007 using exact, simulated data D 1 f. (a),(b) The planes {x 3 = 0.5} and {x 3 = 0.75}. (c),(d) The planes {x 3 = 0.9} and {x 3 = 0.95}.

The very simple vector field f(x) = (1, 0, 0)⊤ allows us to gain more insight into possible problems and limitations of our algorithm. Figure 8 depicts the original vector field and the reconstruction in the plane {x 3 = 0} just as for the other vector fields before. The regularization parameter was γ = 0.0075. There is approximately the same directional error at the boundary as in Figure 4. In addition to that, a slight error in intensity becomes visible. As before, the reconstruction for different planes in the x 3-direction is shown in Figure 9. As in the reconstruction of the field f(x) = (1 − x 2 2 − x 3 2, 0, 0)⊤, it is clearly visible that the directional error at the boundary of the field is reduced and that we obtain a uniform direction of the arrows the farther we move away from the plane {x 3 = 0}. But the images also show that the intensity of the vector field is slightly decreasing with increasing x 3 which should certainly not be the case for this particular vector field. In our ongoing studies of the reconstruction algorithm we try to avoid this problem by either using a varying scaling factor for the different planes or by using a different regularization parameter.

Figure 8

(a) Original vector field f(x) = (1, 0, 0)⊤. (b) Reconstruction with the described algorithm for γ = 0.0075 using exact, simulated data D 1 f.

Figure 9

Reconstruction of the vector field f(x) = (1, 0, 0)⊤ with the described algorithm for γ = 0.0075 using exact, simulated data D 1 f. (a),(b) The planes {x 3 = 0.5} and {x 3 = 0.75}. (c),(d) The planes {x 3 = 0.9} and {x 3 = 0.95}.

Since the described algorithm allows the reconstruction of any arbitrary plane in the x 3-direction, a vertical cross-section of a vector field can easily be computed. Figure 10 shows a reconstruction of the plane {x 1 = 0} for the circular vector field f(x) = (−x 2, x 1, 0)⊤ for two different viewing angles. Figures 11 and 12 display the analogous results for the vector fields f(x) = (1 − x 2 2 − x 3 2, 0, 0)⊤ and f(x) = (1, 0, 0)⊤, respectively. In Figure 11 we recognize a laminar flow just as it would be expected for the flow in blood vessels whereas Figure 12 once more reveals the undesired decrease in intensity.

Figure 10

Reconstruction of the vector field f(x) = (−x 2, x 1, 0)⊤ with the described algorithm for γ = 0.00692 in the plane {x 1 = 0}.

Figure 11

Reconstruction of the vector field f(x) = (1 − x 2 2 − x 3 2, 0, 0)⊤ with the described algorithm for γ = 0.007 in the plane {x 1 = 0}.

Figure 12

Reconstruction of the vector field f(x) = (1, 0, 0)⊤ with the described algorithm for γ = 0.0075 in the plane {x 1 = 0}.

Up to now, we have only shown reconstructions of divergence-free vector fields which were perpendicular to (0, 0, 1)⊤, planar solenoidal fields so to speak. This was done because reconstructing the plane {x 3 = 0} we have to face another difficult problem. Considering the plane {x 3 = 0} it is obvious that the third component of the direction vector ω must always be zero, that is, ω = (ω 1, ω 2, 0)⊤. From formula (3) we can easily calculate the projection of the vector field f onto ω ∈ {x 3 = 0} as which means that the vertical component f 3(x) of any vector field cannot be reconstructed in that particular case. But even for planes where x 3 ≠ 0, the reconstruction of the third component of a vector field will be very difficult. This becomes clear if we look at each of the values ω 1, ω 2 and ω 3 of our direction vector ω = (ω 1, ω 2, ω 3)⊤ separately. Clearly, ω 1 and ω 2 take all values from −1 to 1 if the source travels around the circular scanning curve. This does not apply to ω 3. The maximal value for ω 3 is obtained for the tangential ray with ω 2 = 0, that is the maximal ray in the x 1-x 3-plane hitting the object Ω in just one point. From Figure 13 we easily see that Using the equivalence we obtain For the direction vector ω we can then deduce

Figure 13

Clarifying the calculation of the maximal value for ω 3.

This simple geometric calculation shows that the angle β between the plane {x 3 = 0} and the tangential ray is not exceeding ≈ 19.47° for our geometric setup in which the object is contained in the three-dimensional unit ball Ω3, so R = 1, and the scanning circle Γ has radius r = 3. We conclude that The problem with the small values for ω 3 gets even more difficult the larger the distance between object and scanning circle is chosen. This fact might even be responsible for the intensity error observed when reconstructing different planes in the x 3-direction as mentioned above.

The consequences of the problem can easily be illustrated by applying our algorithm to the two fully three-dimensional vector fields f(x) = (1, 1, 1)⊤ and f(x) = (1, 1, −1)⊤ whose x 3-component should obviously point to opposite directions. Figure 14 shows that even for different perspectives there is no visible difference between the two reconstructions in the plane {x 3 = 0}. Only for larger (or smaller) values of x 3, that means planes above (or below) the plane of the scanning curve {x 3 = 0}, the reconstructions become distinguishable. This is shown in Figure 15 where the results for {x 3 = 0.75} and {x 3 = 0.95} can be directly compared to each other. Even in this case the difference is marginal. For the vector field f(x) = (1, 1, 1)⊤ the arrows characterizing the vectors are slightly pointing upwards whereas those of f(x) = (1, 1, −1)⊤ are pointing downwards. This can be seen from the scale as well as from the colored points indicating the reconstruction plane. For the latter vector field these points are painted above the arrows proving that the reconstructed vectors point in a negative x 3-direction.

Figure 14

(a),(b) Reconstruction of the vector field f(x) = (1, 1, 1)⊤ with the described algorithm for γ = 0.007 in the plane {x 3 = 0}. (c),(d) Reconstruction of the vector field f(x) = (1, 1, −1)⊤ with the described algorithm for γ = 0.007 in the plane {x 3 = 0}.

Figure 15

Top: reconstruction of the vector field f(x) = (1, 1, 1)⊤ with the described algorithm for γ = 0.007 in the planes (a) {x 3 = 0.75} and (b) {x 3 = 0.95}. Bottom: reconstruction of the vector field f(x) = (1, 1, −1)⊤ with the described algorithm for γ = 0.007 in the planes (c) {x 3 = 0.75} and (d) {x 3 = 0.95}.

Despite all these drawbacks it is nevertheless possible to reconstruct a whole three-dimensional vector field. Summarizing the results so far, the x 1- and x 2-component of any vector field f(x) = (f 1, f 2, f 3)⊤, that is, f 1 and f 2, can be reconstructed for any plane in the x 3-direction. Changing our measurement setup by adding to the current scanning curve Γ = r S 2 ∩ {x 3 = 0}, r = 3, a second, orthogonal circle of the same radius in the plane {x 1 = 0} enables us to reconstruct the x 2- and x 3-component of any vector field, that is f 2 and f 3, by using our algorithm as usual. Lade et al. used a comparable setup of two “perpendicular tilt series” in [21] for their longitudinal and transverse measurements. For this additional circle we compute the vertical cross-section of the field in the appropriate plane as shown before. The computations can be done simultaneously which means that no additional time is needed. The modified measurement setup is depicted in Figure 16. As we have shown, we are able to reconstruct f 1 and f 2 for any plane in the x 3-direction with the algorithm presented at the beginning of the paper. The second, orthogonal circle in the plane {x 1 = 0} enables us to analogously reconstruct f 2 and f 3 for any plane in the x 1-direction. Taking a vertical cross-section at x 3 = 0, both results can be combined to reconstruct a complete slice of a three-dimensional vector field. Thereby, we have to pay attention to compute the arithmetic mean for the f 2-component since it is correctly reconstructed for both measurements. First numerical results are depicted in Figure 17, where the vector fields f(x) = (1, 1, 1)⊤ and f(x) = (1, 1, −1)⊤ have been reconstructed by means of the described method. In contrast to the previous images the width of the arrows has been changed to improve the recognizability of the various details.

Figure 16

Measurement setup using two orthogonal circles as scanning curve Γ, one in the plane {x 3 = 0} as usual, and one in the plane {x 1 = 0}.

Figure 17

(a) Reconstruction of the vector field f(x) = (1, 1, 1)⊤ with the described algorithm for γ = 0.007 with circles in the planes {x 3 = 0} and {x 1 = 0}. (b) Reconstruction of the vector field f(x) = (1, 1, −1)⊤ with the described algorithm for γ = 0.007 with circles in the planes {x 3 = 0} and {x 1 = 0}.

It should be added that it is also possible to choose the second orthogonal circle to lie in the plane {x 2 = 0} instead of {x 1 = 0}. This leads to some minor changes in the appearance of the reconstructed images which can be seen from Figure 18. It seems as if each of the vertical circles has its advantages in the reconstruction of a certain direction and as such one or the other may be better suited if some prior knowledge of the vector field exists.

Figure 18

(a) Reconstruction of the vector field f(x) = (1, 1, 1)⊤ with the described algorithm for γ = 0.007 with circles in the planes {x 3 = 0} and {x 2 = 0}. (b) Reconstruction of the vector field f(x) = (1, 1, −1)⊤ with the described algorithm for γ = 0.007 with circles in the planes {x 3 = 0} and {x 2 = 0}.

Moreover, we can combine the advantages of both scanning geometries suited for fully three-dimensional reconstruction we introduced so far. This can be done by simply using data from all three orthogonal circles at once. This obviously leads to a certain redundancy in the data which might nevertheless be useful to improve the quality of our reconstructions. The yellow circle in Figure 19 is not necessary to obtain a complete 3D reconstruction of the vector field, it is only meant to supply us with additional information. Reconstructions of the two vector fields f(x) = (1, 1, 1)⊤ and f(x) = (1, 1, −1)⊤ using the three orthogonal circles can be found in Figure 20. Comparing them with the images made by using only two circles as scanning curve we can see that especially the direction but also the length of the arrows is reconstructed much better.

Figure 19

Measurement setup using all three orthogonal circles as scanning curve Γ, one in each plane {x = 0}, i = 1, 2, 3.

Figure 20

(a) Reconstruction of the vector field f(x) = (1, 1, 1)⊤ with the described algorithm for γ = 0.007 using three orthogonal circles in the planes {x = 0}, i = 1, 2, 3. (b) Reconstruction of the vector field f(x) = (1, 1, −1)⊤ with the described algorithm for γ = 0.007 using three orthogonal circles in the planes {x = 0}, i = 1, 2, 3.

To verify our assumptions we reconstructed once again the planar circular vector field f(x) = (−x 2, x 1, 0)⊤. Comparing the images in Figure 21, where the width of the arrows was reduced in comparison to the ones from the beginning of the paper, we see that using all three orthogonal circles is by far better than using only two of them. Unfortunately we have to admit that the three-dimensional reconstruction is not as good as the one obtained by only one circle in the plane {x 3 = 0}. This can be explained by the fact that the regularization parameter γ was optimized for the latter scanning geometry and was then used for all further modifications of the measurement setup as well. Thus, further enhancements of the results are to be expected by choosing a varying regularization parameter for the different circular scanning geometries. In addition to that we might implement a scaling factor to compensate for the loss of intensity at the boundary of our reconstruction region. It might also be possible to incorporate some sort of prior knowledge. For example in clinical applications when measuring blood flow we may use some information about the blood vessels to correct the direction of the flow at the boundary.

Figure 21

(a) Original vector field f(x) = (−x 2, x 1, 0)⊤. (b) Reconstruction with the described algorithm for γ = 0.00692 using one single circle in the plane {x 3 = 0}. (c) Reconstruction with the described algorithm for γ = 0.00692 using circles in the planes {x 3 = 0} and {x 1 = 0}. (d) Reconstruction with the described algorithm for γ = 0.00692 using circles in the planes {x 3 = 0} and {x 2 = 0}. (e) Reconstruction with the described algorithm for γ = 0.00692 using three orthogonal circles in the planes {x = 0}, i = 1, 2, 3.

It is to mention that in contrast to our initial scanning geometry of only one circle the two as well as the three orthogonal circles certainly meet the requirements of Tuy-Kirillov's condition [27, 28], namely that the source curve Γ intersects each plane hitting supp(f) transversally, see [17]. This might be useful for future theoretical advances in the field of three-dimensional vector tomography as well as for improvements and extensions to our algorithm. This is part of our current research.

4. CONCLUSION

We presented a first approach for reconstructing cone beam data in three-dimensional vector field tomography. The algorithm relies on known results for the scalar case from [12, 26], where the method of approximate inverse has been applied for the computation of reconstruction kernels for circular 3D cone beam tomography. A possibility to extend that very efficient regularization technique to VFT has been shown and first numerical experiments are very promising.

The investigation of what happens when we use a scanning curve Γ different from what has been presented in this paper is subject of current research. Nevertheless, it is to mention that the proposed algorithm is not restricted to circles as scanning curves. Helical scanning geometries are especially interesting in clinical applications because a helical source trajectory can easily be implemented by moving the patient's bed through the scanner's gantry at constant speed, a method which is common practice in cone beam CT scanning today. Although data acquisition itself is more difficult, it is much faster and thus better suited for clinical diagnostics. The results of Katsevich [29] could help when we consider a helix as trajectory. Reconsidering the generalizations to tensor fields of arbitrary rank and in arbitrary dimension from which we refrained in the first section of the paper, it is to mention that Denisjuk [30] formulated a generalization of Tuy-Kirillov's condition to tensor fields of rank m which might help to compute exact reconstruction kernels for D . Denisjuk also proved in [30] that a full reconstruction of solenoidal vector fields is possible if the vectorial cone beam data D 1 f(α, ω) are available for all α on a trajectory Γ satisfying a generalized Tuy condition and all directions ω ∈ S 2.

Furthermore, the problem of signal attenuation which is important for Doppler tomography using ultrasound has been addressed by several authors. Unfortunately, so far only works have been published dealing with the problem in two dimensions. Bukhgeim and Kazantsev derived a 2D inversion formula in [31]. Their proof uses a coordinate transformation into complex variables which cannot be generalized to three dimensions. In [32], Natterer develops an extension of Novikov's inversion formula for 2D vector fields. Finally, Stråhlén proves a Fourier slice theorem for the attenuated vectorial Radon transform in two dimensions for the parallel geometry in [33]. Further advances in this field of research will certainly be interesting for 3D cone beam vector field tomography as well.

By now the data acquisition in vector field tomography neglects the wave structure of the ultrasound signals. Further research might take the wave structure into account to improve the reconstruction and modeling accuracy. In this respect the application of the Rayleigh-Sommerfeld formula could be useful.