A Reconstruction Approach for Imaging in 3D Cone Beam Vector Field Tomography

Abstract

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 {𝑥₃=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 {𝑥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 𝑚 by𝐃_𝑚𝐟(𝛼,𝜔)=^∞₀⟨𝐟(𝛼+𝑡𝜔),𝜔^𝑚=⟩d𝑡_∞0𝐟_{𝑖1⋯𝑖𝑚}(𝛼+𝑡𝜔)𝜔^𝑖1⋯𝜔^𝑖𝑚d𝑡, (1) where 𝛼∈Γ is a source point on the scanning curve Γ⊂ℝ^𝑛⧵Ω which surrounds the object Ω,𝜔∈𝑆^𝑛−1 is the unit vector of direction of the line and 𝐟 is a tensor field of rank 𝑚 with compact support in the open domain Ω. Tensor fields of rank 𝑚=0 are scalar functions 𝑓(𝑥), tensor fields of rank 𝑚=1 are vector fields 𝐟(𝑥) in ℝ^𝑛. In (1) we use Einstein’s summation rule, that means we sum up over equal indices 𝑖_𝑗, where 1≤𝑖_𝑗≤𝑛.

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 𝑛 dimensions. In the following, scalar fields will be denoted by 𝑓 whereas vector fields will be written bold-faced, such as 𝐟.

Setting 𝑚=0 in (1) we obtain the well-known (scalar) cone beam transform𝐃₀𝑓(𝛼,𝜔)=^∞₀𝑓(𝛼+𝑡𝜔)d𝑡 (2)of a scalar field 𝑓∶ℝ^𝑛⊃Ω→ℝ. For 𝑚=1 formula (1) is the cone beam transform for vector fields𝐟∶ℝ^𝑛⊃Ω→ℝ^𝑛, which in 3D is also often referred to as the Doppler transform. It reads as 𝐃₁𝐟(𝛼,𝜔)=^∞₀⟨𝐟(𝛼+𝑡𝜔),𝜔⟩d𝑡. (3) Hence the mathematical problem of 3D cone beam VFT consists of inverting 𝐃₁𝐟=𝑔 for given measurements 𝑔∈ℝ, that is reconstructing a three-dimensional vector field 𝐟 from one-dimensional, that is, scalar, data 𝑔. In contrast to 𝐃₀ an inversion formula for 𝐃₁ 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 𝐟(𝑥)=∇𝜙(𝑥)=(𝜕_𝑥1𝜙(𝑥),…,𝜕_𝑥𝑛𝜙(𝑥))^⊤,𝜙∶ℝ^𝑛⊃Ω→ℝ, with compact support 𝐃₁𝐟(𝛼,𝜔)=𝐃₁(∇𝜙)(𝛼,𝜔)=^∞₀∇𝜙(𝛼+𝑡𝜔)𝜔d𝑡=[𝜙(𝛼+𝑡𝜔)]^∞₀=0. (4)

Figure 2 

Visualization of the measured projection data for 𝐃1𝐟.

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 𝐟 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 𝑔 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 𝑚=0 and 𝑚=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 𝑛=3,𝑚=0 (see (2)) in [14]. Then we outline how the method of approximate inverse can be used for (2) to solve 𝐃₀𝑓=𝑔 and present its application to (3) to solve 𝐃₁𝐟=𝑔. Furthermore, an approach is presented, how reconstruction kernels for 𝐃₁ from (3) can be calculated with the help of known reconstruction kernels for the scalar cone beam transform 𝐃₀. Some pictures from numerical experiments show that this approach is very promising.

2. Mathematical Properties of 𝐃₀ and 𝐃₁

Let 𝐿²(𝑋,ℝ^𝑛)∶={𝐟∶𝑋→ℝ^𝑛∶‖𝐟‖_𝐿²=⟨𝐟,𝐟⟩_𝐿^1/2₂<∞} (5) denote the space of square integrable functions from 𝑋⊂ℝ^𝑛 to ℝ^𝑛 where the 𝐿²-inner product of two functions is given as ⟨𝐟,𝐠⟩_𝐿²=_𝑋⟨𝐟(𝑥),𝐠(𝑥)⟩_ℝ^𝑛d𝑥. (6) Fundamental properties of (2) and (3) are summarized in the following theorem, which is the result of straightforward calculations.

Theorem 1. Let Ω^𝑛∶={𝑥∈ℝ^𝑛∶|𝑥|<1} with 𝜕Ω^𝑛=𝑆^𝑛−1 be the reconstruction region, that is, the region in which the object Ω is contained. The mappings𝐃₀∶𝐿²(Ω^𝑛)→𝐿²(Γ×𝑆^𝑛−1𝐃),₁∶𝐿²(Ω^𝑛,ℝ^𝑛)→𝐿²(Γ×𝑆^𝑛−1) (7) are linear and bounded if _Γ(|𝛼|−1)^1−𝑛𝑑𝛼<∞. (8)The adjoints (backprojections) 𝐃^∗₀∶𝐿²(Γ×𝑆^𝑛−1)→𝐿²(Ω^𝑛) and 𝐃^∗₁∶𝐿²(Γ×𝑆^𝑛−1)→𝐿²(Ω^𝑛,ℝ^𝑛) are given by 𝐃^∗₀𝑔(𝑥)=_Γ{|𝑥−𝛼|^1−𝑛𝑔(𝛼,𝑥−𝛼𝐃|𝑥−𝛼|)}𝑑𝛼,_∗1𝑔(𝑥)=_Γ{|𝑥−𝛼|^−𝑛𝑔(𝛼,𝑥−𝛼|𝑥−𝛼|)(𝑥−𝛼)}𝑑𝛼. (9)

For 𝑛=3, we obtain the well-known cone beam transform 𝐃₀ with the corresponding backprojection operator 𝐃^∗₀ 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𝛼. (10)𝐃^∗₁𝑔(𝑥) represents an integration over all lines intersecting 𝑥. The result can be seen as an average over all directions connecting a source point 𝛼 and 𝑥.

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 𝑚 in 𝑛 dimensions in [13] but our presentation will be restricted to 𝐃₀ and 𝐃₁.

Theorem 2 (Schuster [13] based on Hamaker et al. [15]). Assume 𝑛≥2,𝑓∈𝒞₀^(𝑛−2)(Ω^𝑛), and 𝐟∈𝒞₀^(𝑛−2)(Ω^𝑛,ℝ^𝑛). Then, 𝜕^(𝑛−2)𝜕𝑠^(𝑛−2)𝐑𝑓(𝜔,⟨𝛼,𝜔⟩)=(−1)^(𝑛−2)_𝑆𝑛−1𝐃₀𝑓(𝛼,𝜃)𝛿^(𝑛−2)𝜕(⟨𝜔,𝜃⟩)𝑑𝑆(𝜃),_(𝑛−2)𝜕𝑠^(𝑛−2)𝐑𝑓_𝛼(𝜔,⟨𝛼,𝜔⟩)=(−1)^(𝑛−2)_𝑆𝑛−1𝐃₁𝐟(𝛼,𝜃)𝛿^(𝑛−2)(⟨𝜔,𝜃⟩)𝑑𝑆(𝜃), (11)(12) where 𝛼∈Γ,𝜔∈𝑆^𝑛−1,𝑑𝑆 denotes the surface measure on 𝑆^𝑛−1,𝐑 is the 𝑛-dimensional Radon transform defined by(𝐑𝑓)(𝜃,𝑠)=_{⟨𝑥,𝜃⟩=𝑠}𝑓(𝑥)𝑑𝑥, (13)and𝑓_𝛼(𝑥)=⟨𝐟(𝑥),|𝑥−𝛼|⁻¹(𝑥−𝛼)⟩_ℝ^𝑛 (14)is the projection of 𝐟 onto (𝑥−𝛼)/|𝑥−𝛼|.

In the following, we want to limit our presentation to 𝑛=3. Then formula (12) reads as 𝜕𝜕𝑠𝐑𝑓_𝛼(𝜔,⟨𝛼,𝜔⟩)=(−1)_𝑆2𝐃₁𝐟(𝛼,𝜃)𝛿^′=(⟨𝜔,𝜃⟩)d𝑆(𝜃)_{𝑆2∩{⟨𝜃,𝜔⟩=0}}⟨∇_𝑦𝐃₁𝐟(𝛼,𝑦=𝜃),𝜔⟩d𝑆(𝜃). (15)

In the scalar case a solver for 𝐃₀ can be constructed with the help of (11) for 𝑛=3. This is done by inverting the Radon transform 𝐑 which is possible if the condition of Tuy-Kirillov is satisfied. It tells that we have full knowledge of 𝐑𝑔(𝜔,𝑠) for all 𝜔,𝑠 and any scalar function 𝑔∶Ω³→ℝ, if any plane intersecting the object Ω⊆Ω³ does also have at least one intersection point with the scanning curve Γ and this intersection must be nontransversally. This works fine for 𝐃₀ since then 𝑓(𝑥) is independent of 𝛼 but unfortunately that does not help in case of 𝐃₁ (and analogous transforms for tensor fields as well, see [13]), since there the object function 𝑓_𝛼=⟨𝐟(𝑥),(𝑥−𝛼)/|𝑥−𝛼|⟩_ℝ³ of 𝐑 depends on and hence changes with 𝛼, see (15). Thus we seek an alternative way of solving 𝐃₁𝐟=𝑔 for vector fields 𝐟.

3. Approximation of Reconstruction Kernels in Vector Field Tomography

As already said, an inversion formula for 𝐃₁ 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 𝐃₁. 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 𝐃₁. 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 𝐃₀. We describe this approach and then formulate an extension of it to 𝐃₁.

Let 𝑓∈𝐿²(Ω³) be a scalar function. The approximate inverse computes a smoothed version 𝑓_𝛾 of 𝑓 by convolving 𝑓 with a mollifier𝑒_𝛾∈𝒞^∞(ℝ³). A mollifier 𝑒_𝛾 is a smooth function with small essential support having the property that 𝑓_𝛾(𝑥)∶=(𝑓∗𝑒_𝛾)(𝑥)→𝑓(𝑥)as𝛾→0. (16)Here, ∗ denotes the convolution(𝑓∗ℎ)(𝑥)=_ℝ3𝑓(𝑦−𝑥)ℎ(𝑦)d𝑦. (17) Such a function is given by the Gaussian kernel 𝑒_𝛾𝛾(𝑥)=⁻³(2𝜋)^3/2exp(−|𝑥|²2𝛾²). (18) Provided that we can solve the equation 𝐃^∗₀[𝑣_𝛾(𝑥)]=𝑒_𝛾(𝑥−⋅), (19) then we can reconstruct 𝑓_𝛾 from the measured cone beam data 𝐃₀𝑓 by 𝑓_𝛾(𝑥)=⟨𝐃₀𝑓,𝑣_𝛾(𝑥)⟩_{𝐿2(Γ×𝑆2)}=_Γ𝑆2(𝐃₀𝑓)(𝛼,𝜔)𝑣_𝛾(𝑥;𝛼,𝜔)d𝑆(𝜔)d𝛼, (20) where 𝑣_𝛾(𝑥)=𝑣_𝛾(𝑥;𝛼,𝜔)∈𝐿²(Γ×𝑆²) for 𝑥∈Ω^𝑛 is called a reconstruction kernel. Hence the method of approximate inverse consists of evaluating inner products of the given data 𝐃₀𝑓 with reconstruction kernels 𝑣_𝛾(𝑥). This can be done in an efficient way using the translation invariance of 𝑒_𝛾 and some appropriate approximations which are outlined in detail in [12]. These imply that we have to solve (19) only once, namely for 𝑥=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 𝐃₁ and hence to VFT, we construct mollifier fields 𝐄^𝑗_𝛾∈𝐿²(Ω³,ℝ³) defining 𝐄^𝑗_𝛾(𝑥)∶=𝑒_𝛾(𝑥)⋅𝑒_𝑗,𝑗∈{1,2,3}, (21) where 𝑒₁=(1,0,0)^⊤,𝑒₂=(0,1,0)^⊤ and 𝑒₃=(0,0,1)^⊤. Using again the Gaussian (18) as mollifier 𝑒_𝛾 we obtain (𝐟_𝛾)_𝑗(𝑥)∶=(𝐟∗𝐄^𝑗_𝛾)(𝑥)→𝐟_𝑗(𝑥)as𝛾→0 (22) for 𝐟∈𝐿²(Ω³,ℝ³). Unfortunately, by now the exact reconstruction kernels 𝐕^𝑗_𝛾(𝑥), that is, the solutions of 𝐃^∗₁[𝐕^𝑗_𝛾(𝑥)]=𝐄^𝑗_𝛾(𝑥−⋅) are still unknown. But the special structure of the mollifier fields 𝐄^𝑗_𝛾 allow for a computation of reconstruction kernels for 𝐏𝐟(𝛼,𝜔)=^∞₀𝐟(𝛼+𝑡𝜔)d𝑡 (23) with the help of kernels for 𝐃₀.

Theorem 3. Let 𝑣_𝛾 be the reconstruction kernel associated to 𝑒_𝛾 with respect to 𝐃₀, that is, 𝐃^∗₀[𝑣_𝛾(𝑥)]=𝑒_𝛾(𝑥−⋅). (24) Defining 𝐕^𝑗_𝛾(𝑥;𝛼,𝜔)∶=𝑣_𝛾(𝑥;𝛼,𝜔)⋅𝑒_𝑗∈𝐿²(Γ×𝑆²,ℝ³) yields 𝐏^∗[𝐕^𝑗_𝛾(𝑥)]=𝐄^𝑗_𝛾(𝑥−⋅), (25) that means 𝐕^𝑗_𝛾 is a reconstruction kernel associated to 𝐄^𝑗_𝛾 with respect to 𝐏. The adjoint 𝐏^∗ of 𝐏 is given as 𝐏^∗𝐠(𝑥)=_Γ|𝑥−𝛼|⁻²𝐠(𝛼,𝑥−𝛼|𝑥−𝛼|)𝑑𝛼 (26) for 𝐠∈𝐿²(Γ×𝑆²,ℝ³).

Proof. The adjoint 𝐏^∗ is computed as ⟨𝐏𝐟,𝐠⟩_𝐿²_(Γ×𝑆²_,ℝ³₎=_Γ𝑆2^∞₀=⟨𝐟(𝛼+𝑡𝜔),𝐠(𝛼,𝜔)⟩d𝑡d𝑆(𝜔)d𝛼_Γ𝑅3⟨𝐟(𝑥),𝐠(𝛼,𝑥−𝛼|𝑥−𝛼|)⟩|𝑥−𝛼|⁻²d𝑥d𝛼=⟨𝐟,𝐏^∗𝐠⟩_𝐿²_(Ω,ℝ³₎, (27) where we applied Fubini’s theorem as well as the substitution 𝑥=𝛼+𝑡𝜔. A short calculation further shows that𝐏^∗[𝐕^𝑗_𝛾(𝑥)](𝑦)=_Γ|𝑦−𝛼|⁻²𝐕^𝑗_𝛾(𝑥;𝛼,𝑦−𝛼=|𝑦−𝛼|)d𝛼_Γ|𝑦−𝛼|⁻²𝑣_𝛾(𝑥;𝛼,𝑦−𝛼|𝑦−𝛼|)d𝛼⋅𝑒_𝑗=𝐃^∗₀[𝑣_𝛾(𝑥)](𝑦)⋅𝑒_𝑗=𝑒_𝛾(𝑥−𝑦)⋅𝑒_𝑗=𝐄^𝑗_𝛾(𝑥−𝑦). (28)

The data 𝐏𝐟 are not known and cannot be computed from 𝐃₁𝐟. But, observing that 𝐃₁𝐟(𝛼,𝜔)=⟨𝐏𝐟(𝛼,𝜔),𝜔⟩_ℝ³ (29) and since 𝐏𝐟(𝛼,𝜔)∈ℝ³ we obtain 𝐏𝐟(𝛼,𝜔)=𝐃₁𝐟(𝛼,𝜔)𝜔+𝜆₁(𝛼,𝜔^⟂₁)𝜔^⟂₁+𝜆₂(𝛼,𝜔^⟂₂)𝜔^⟂₂, (30) where 𝜔^⟂₁,𝜔^⟂₂∈𝑆² are such that {𝜔,𝜔^⟂₁,𝜔^⟂₂} is an orthonormal basis of ℝ³ and 𝜆₁,𝜆₂ are appropriate coefficients. Thus approximating 𝐏𝐟(𝛼,𝜔)≈𝐃₁𝐟(𝛼,𝜔)𝜔 (31) we neglect the parts orthogonal to 𝜔 and can apply the method of approximate inverse using the reconstruction kernels 𝐕^𝑗_𝛾 for 𝐏. This procedure results in Algorithm 1.

Given   : Measured data 𝐃1𝐟(𝛼,𝜔) for 𝛼∈Γ, 𝜔∈𝑆2 Output    : Approximation 𝐟𝛾 to 𝐟 Compute:     (1) 𝐠(𝛼,𝜔)=𝐃1𝐟(𝛼,𝜔)𝜔     (2) (𝐟𝛾)𝑗(𝑥)=⟨𝐠,𝐕𝑗𝛾(𝑥)⟩𝐿2(Γ×𝑆2,ℝ3)            =∫Γ∫𝑆2⟨𝐠(𝛼,𝜔),𝐕𝑗𝛾(𝑥;𝛼,𝜔)⟩d𝑆(𝜔)d𝛼    =∫Γ∫𝑆2𝐃1𝐟(𝛼,𝜔)𝜔𝑗𝑣𝛾(𝑥;𝛼,𝜔)d𝑆(𝜔)d𝛼                       for 𝑗∈{1,2,3}

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 𝐃₀𝑓(𝛼,𝜔) are transferred to its 3D equivalent 𝐏𝐟(𝛼,𝜔) which is an approximation to 𝐃₁𝐟(𝛼,𝜔) by using 𝐏𝐟(𝛼,𝜔)≈𝐃₁𝐟(𝛼,𝜔)𝜔.

Figures 3 and 4 display first results of the above algorithm when applied to exact, simulated data for the solenoidal vector fields 𝐟(𝑥)=(−𝑥₂,𝑥₁,0)^⊤ and 𝐟(𝑥)=(1−𝑥²₂−𝑥²₃,0,0)^⊤. As already said, we made use of the measurement setup as shown in Figure 1. The scanning curve was Γ=𝑟𝑆²∩{𝑥₃=0},𝑟=3, that is a circle of radius 𝑟=3 in the plane {𝑥₃=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 Ω=Ω³. The reconstructions depicted in Figures 3 and 4 were done in the plane {𝑥₃=0}. The mollifier 𝑒_𝛾 defining the fields 𝐄^𝑗_𝛾 was chosen as the Gaussian given in (18). The regularization parameter was 𝛾=0.00692 and 𝛾=0.007, respectively. The corresponding reconstruction kernel 𝑣_𝛾(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 𝐟 itself. Nevertheless, in our experiments a value of 𝛾≈0.007 always led to good results.

(a)

(b)

Figure 3 

(a) Original vector field 𝐟(𝑥)=(−𝑥2,𝑥1,0)⊤. (b) Reconstruction with the described algorithm for 𝛾=0.00692 using exact, simulated data 𝐃1𝐟.

(a)

(b)

Figure 4 

(a) Original vector field 𝐟(𝑥)=(1−𝑥22−𝑥23,0,0)⊤. (b) Reconstruction with the described algorithm for 𝛾=0.007 using exact, simulated data 𝐃1𝐟.

Figure 5 

Reconstruction kernel 𝑣𝛾(0;𝛼,𝜔) for 𝛾=0.007 associated to the scalar cone beam transform 𝐃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 𝑥₃-direction. Figures 6 and 7 show the reconstructions for the planes {𝑥₃=0.5},{𝑥₃=0.75},{𝑥₃=0.9} and {𝑥₃=0.95} of the two afore-mentioned vector fields. Figure 7 also illustrates that the intensity of the vector field 𝐟(𝑥)=(1−𝑥²₂−𝑥²₃,0,0)^⊤ decreases as should be expected since 𝑥²₃ is subtracted in the first component and that even the directional error at the boundary is reduced with increasing 𝑥₃.

(a)

(b)

(c)

(d)

Figure 6 

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

(a)

(b)

(c)

(d)

Figure 7 

Reconstruction of the vector field 𝐟(𝑥)=(1−𝑥22−𝑥23,0,0)⊤ with the described algorithm for 𝛾=0.007 using exact, simulated data 𝐃1𝐟. (a),(b) The planes {𝑥3=0.5} and {𝑥3=0.75}. (c),(d) The planes {𝑥3=0.9} and {𝑥3=0.95}.

The very simple vector field 𝐟(𝑥)=(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 {𝑥₃=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 𝑥₃-direction is shown in Figure 9. As in the reconstruction of the field 𝐟(𝑥)=(1−𝑥²₂−𝑥²₃,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 {𝑥₃=0}. But the images also show that the intensity of the vector field is slightly decreasing with increasing 𝑥₃ 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.

(a)

(b)

Figure 8 

(a) Original vector field 𝐟(𝑥)=(1,0,0)⊤. (b) Reconstruction with the described algorithm for 𝛾=0.0075 using exact, simulated data 𝐃1𝐟.

(a)

(b)

(c)

(d)

Figure 9 

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

Since the described algorithm allows the reconstruction of any arbitrary plane in the 𝑥₃-direction, a vertical cross-section of a vector field can easily be computed. Figure 10 shows a reconstruction of the plane {𝑥₁=0} for the circular vector field 𝐟(𝑥)=(−𝑥₂,𝑥₁,0)^⊤ for two different viewing angles. Figures 11 and 12 display the analogous results for the vector fields 𝐟(𝑥)=(1−𝑥²₂−𝑥²₃,0,0)^⊤ and 𝐟(𝑥)=(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.

(a)

(b)

Figure 10 

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

(a)

(b)

Figure 11 

Reconstruction of the vector field 𝐟(𝑥)=(1−𝑥22−𝑥23,0,0)⊤ with the described algorithm for 𝛾=0.007 in the plane {𝑥1=0}.

(a)

(b)

Figure 12 

Reconstruction of the vector field 𝐟(𝑥)=(1,0,0)⊤ with the described algorithm for 𝛾=0.0075 in the plane {𝑥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 {𝑥₃=0} we have to face another difficult problem. Considering the plane {𝑥₃=0} it is obvious that the third component of the direction vector 𝜔 must always be zero, that is, 𝜔=(𝜔₁,𝜔₂,0)^⊤. From formula (3) we can easily calculate the projection of the vector field 𝐟 onto 𝜔∈{𝑥₃=0} as 𝐃₁𝐟(𝛼,𝜔)=^∞₀=⟨𝐟(𝛼+𝑡𝜔),𝜔⟩d𝑡_∞0(𝑓₁(𝛼+𝑡𝜔)𝜔₁+𝑓₂(𝛼+𝑡𝜔)𝜔₂+𝑓₃(𝛼+𝑡𝜔)𝜔₃)d𝑡_𝜔3=0=^∞₀(𝑓₁(𝛼+𝑡𝜔)𝜔₁+𝑓₂(𝛼+𝑡𝜔)𝜔₂)d𝑡, (32) which means that the vertical component 𝑓₃(𝑥) of any vector field cannot be reconstructed in that particular case. But even for planes where 𝑥₃≠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 𝜔₁,𝜔₂ and 𝜔₃ of our direction vector 𝜔=(𝜔₁,𝜔₂,𝜔₃)^⊤ separately. Clearly, 𝜔₁ and 𝜔₂ take all values from −1 to 1 if the source travels around the circular scanning curve. This does not apply to 𝜔₃. The maximal value for 𝜔₃ is obtained for the tangential ray with 𝜔₂=0, that is the maximal ray in the 𝑥₁-𝑥₃-plane hitting the object Ω in just one point. From Figure 13 we easily see that 𝑅sin(𝛽)=𝑟,𝑢tan(𝛽)=𝑟𝑅⟺𝑢=𝑟tan(𝛽)=𝑟tan(arcsin(𝑟)). (33) Using the equivalence 𝑥arcsin(𝑥)=arctan(√1−𝑥²) (34) we obtain 𝑢=𝑟tan(arctan(𝑅/𝑟1−(𝑅/𝑟)²=𝑅))(𝑟²−𝑅²)/𝑟²=𝑟𝑅𝑟²−𝑅². (35) For the direction vector 𝜔 we can then deduce 1𝜔=0𝑢,√‖𝜃‖𝜃with𝜃=−𝑟‖𝜃‖=𝑟²+𝑢²=𝑟²+(𝑟𝑅)²𝑟²−𝑅²=𝑟𝑅1+²𝑟²−𝑅²=𝑟𝑟²−𝑅²+𝑅²𝑟²−𝑅²=𝑟²𝑟²−𝑅²,𝑢√𝑟²+𝑢²=(𝑟𝑅/𝑟²−𝑅²)(𝑟²/𝑟²−𝑅²)=𝑅𝑟,−𝑟√𝑟²+𝑢²=−𝑟(𝑟²/𝑟²−𝑅²)(=−𝑟²−𝑅²)𝑟−(⇒𝜔=𝑟²−𝑅²)𝑟0𝑅𝑟. (36)

Figure 13 

Clarifying the calculation of the maximal value for 𝜔3.

This simple geometric calculation shows that the angle 𝛽 between the plane {𝑥₃=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 Ω³, so 𝑅=1, and the scanning circle Γ has radius 𝑟=3. We conclude that −𝑅𝑟≤𝜔₃≤𝑅𝑟1⟺−3≤𝜔₃≤13. (37) The problem with the small values for 𝜔₃ 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 𝑥₃-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 𝐟(𝑥)=(1,1,1)^⊤ and 𝐟(𝑥)=(1,1,−1)^⊤ whose 𝑥₃-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 {𝑥₃=0}. Only for larger (or smaller) values of 𝑥₃, that means planes above (or below) the plane of the scanning curve {𝑥₃=0}, the reconstructions become distinguishable. This is shown in Figure 15 where the results for {𝑥₃=0.75} and {𝑥₃=0.95} can be directly compared to each other. Even in this case the difference is marginal. For the vector field 𝐟(𝑥)=(1,1,1)^⊤ the arrows characterizing the vectors are slightly pointing upwards whereas those of 𝐟(𝑥)=(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 𝑥₃-direction.

(a)

(b)

(c)

(d)

Figure 14 

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

(a)

(b)

(c)

(d)

Figure 15 

Top: reconstruction of the vector field 𝑓(𝑥)=(1,1,1)⊤ with the described algorithm for 𝛾=0.007 in the planes (a) {𝑥3=0.75} and (b){𝑥3=0.95}. Bottom: reconstruction of the vector field 𝑓(𝑥)=(1,1,−1)⊤ with the described algorithm for 𝛾=0.007 in the planes (c) {𝑥3=0.75} and (d) {𝑥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 𝑥₁- and 𝑥₂-component of any vector field 𝐟(𝑥)=(𝑓₁,𝑓₂,𝑓₃)^⊤, that is, 𝑓₁ and 𝑓₂, can be reconstructed for any plane in the 𝑥₃-direction. Changing our measurement setup by adding to the current scanning curve Γ=𝑟𝑆²∩{𝑥₃=0},𝑟=3, a second, orthogonal circle of the same radius in the plane {𝑥₁=0} enables us to reconstruct the 𝑥₂- and 𝑥₃-component of any vector field, that is 𝑓₂ and 𝑓₃, 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 𝑓₁ and 𝑓₂ for any plane in the 𝑥₃-direction with the algorithm presented at the beginning of the paper. The second, orthogonal circle in the plane {𝑥₁=0} enables us to analogously reconstruct 𝑓₂ and 𝑓₃ for any plane in the 𝑥₁-direction. Taking a vertical cross-section at 𝑥₃=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 𝑓₂-component since it is correctly reconstructed for both measurements. First numerical results are depicted in Figure 17, where the vector fields 𝐟(𝑥)=(1,1,1)^⊤ and 𝐟(𝑥)=(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 {𝑥3=0} as usual, and one in the plane {𝑥1=0}.

(a)

(b)

Figure 17 

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

It should be added that it is also possible to choose the second orthogonal circle to lie in the plane {𝑥₂=0} instead of {𝑥₁=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.

(a)

(b)

Figure 18 

(a) Reconstruction of the vector field 𝐟(𝑥)=(1,1,1)⊤ with the described algorithm for 𝛾=0.007 with circles in the planes {𝑥3=0} and {𝑥2=0}. (b) Reconstruction of the vector field 𝐟(𝑥)=(1,1,−1)⊤ with the described algorithm for 𝛾=0.007 with circles in the planes {𝑥3=0} and {𝑥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 𝐟(𝑥)=(1,1,1)^⊤ and 𝐟(𝑥)=(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 {𝑥𝑖=0},𝑖=1,2,3.

(a)

(b)

Figure 20 

(a) Reconstruction of the vector field 𝐟(𝑥)=(1,1,1)⊤ with the described algorithm for 𝛾=0.007 using three orthogonal circles in the planes {𝑥𝑖=0}, 𝑖=1,2,3. (b) Reconstruction of the vector field 𝐟(𝑥)=(1,1,−1)⊤ with the described algorithm for 𝛾=0.007 using three orthogonal circles in the planes {𝑥𝑖=0},𝑖=1,2,3.

To verify our assumptions we reconstructed once again the planar circular vector field 𝐟(𝑥)=(−𝑥₂,𝑥₁,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 {𝑥₃=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.

(a)

(b)

(c)

(d)

(e)

Figure 21 

(a) Original vector field 𝐟(𝑥)=(−𝑥2,𝑥1,0)⊤. (b) Reconstruction with the described algorithm for 𝛾=0.00692 using one single circle in the plane {𝑥3=0}. (c) Reconstruction with the described algorithm for 𝛾=0.00692 using circles in the planes {𝑥3=0} and {𝑥1=0}. (d) Reconstruction with the described algorithm for 𝛾=0.00692 using circles in the planes {𝑥3=0} and {𝑥2=0}. (e) Reconstruction with the described algorithm for 𝛾=0.00692 using three orthogonal circles in the planes {𝑥𝑖=0},𝑖=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(𝑓) 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 𝑚 which might help to compute exact reconstruction kernels for 𝐃_𝑚. Denisjuk also proved in [30] that a full reconstruction of solenoidal vector fields is possible if the vectorial cone beam data 𝐃₁𝐟(𝛼,𝜔) are available for all 𝛼 on a trajectory Γ satisfying a generalized Tuy condition and all directions 𝜔∈𝑆².

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.

Acknowledgment

This project is being supported by the Deutsche Forschungsgemeinschaft (DFG) under Grant Schu 1978/1-6.