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

Authors: T. Schuster; D. Theis; A. K. Louis
Journal: International Journal of Biomedical Imaging
 (2008)

---

## 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.


---

## Body

## 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 {x3=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{x3=0} as scanning
curve Γ.As in scalar 3D computerized tomography, thecone
beam transform is of special interest from a practical point of view. It is
defined for a tensor field of rank m by(1)Dmf(α,ω)=∫0∞〈f(α+tω),ωm〉dt=∫0∞fi1⋯im(α+tω)ωi1⋯ωimdt,
where α∈Γ is a source point on the scanning curve Γ⊂ℝn∖Ω¯ which surrounds the object Ω,ω∈Sn−1 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 ℝn. In (1) we use Einstein’s summation rule, that means we sum up over equal
indices ij, where 1≤ij≤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 thescalar 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.Settingm=0 in (1) we obtain the well-known (scalar)
cone beam transform(2)D0f(α,ω)=∫0∞f(α+tω)dtof a scalar field f:ℝn⊃Ω→ℝ. For m=1 formula (1) is the cone beam transform for
vector fieldsf:ℝn⊃Ω→ℝn, which in 3D is also often referred to as the Doppler transform. It reads as
(3)D1f(α,ω)=∫0∞〈f(α+tω),ω〉dt.
Hence the mathematical problem
of 3D cone beam VFT consists of inverting D1f=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 D0 an inversion formula for D1 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)=(∂x1ϕ(x),…,∂xnϕ(x))⊤,ϕ:ℝn⊃Ω→ℝ, with compact support
(4)D1f(α,ω)=D1(∇ϕ)(α,ω)=∫0∞∇ϕ(α+tω)ωdt=[ϕ(α+tω)]0∞=0.Figure 2
Visualization of the measured projection data forD1f.The method ofapproximate 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 D0f=g and present its application to (3) to solve D1f=g. Furthermore, an approach is presented, how reconstruction kernels for D1 from (3) can be calculated with the help of
known reconstruction kernels for the scalar cone beam transform D0. Some pictures from numerical experiments show that this approach is very promising.


## 2. Mathematical Properties ofD0 and D1

Let(5)L2(X,ℝn):={f:X→ℝn:∥f∥L2=〈f,f〉L21/2<∞}
denote the space of square
integrable functions from X⊂ℝn to ℝn where the L2-inner product of two functions is given as
(6)〈f,g〉L2=∫X〈f(x),g(x)〉ℝndx.
Fundamental properties of (2) and (3) are summarized in the following theorem, which is the result of
straightforward calculations.Theorem 1.
LetΩn:={x∈ℝn:|x|<1} with ∂Ωn=Sn−1 be the reconstruction region, that is, the
region in which the object Ω is contained. The mappings(7)D0:L2(Ωn)→L2(Γ×Sn−1),D1:L2(Ωn,ℝn)→L2(Γ×Sn−1)
are linear and bounded if
(8)∫Γ(|α|−1)1−ndα<∞.The adjoints (backprojections) D0∗:L2(Γ×Sn−1)→L2(Ωn) and D1∗:L2(Γ×Sn−1)→L2(Ωn,ℝn) are given by
(9)D0∗g(x)=∫Γ{|x−α|1−ng(α,x−α|x−α|)}dα,D1∗g(x)=∫Γ{|x−α|−ng(α,x−α|x−α|)(x−α)}dα.Forn=3, we obtain the well-known cone beam transform D0 with the corresponding backprojection operator D0∗ 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
(10)D1∗g(x)=∫Γ|x−α|−2g(α,x−α|x−α|)x−α|x−α|dα.D1∗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 D0 and D1.Theorem 2  (Schuster [13] based on Hamaker et al. [15]).
Assumen≥2,f∈𝒞0(n−2)(Ωn), and f∈𝒞0(n−2)(Ωn,ℝn). Then,
(11)∂(n−2)∂s(n−2)Rf(ω,〈α,ω〉)=(−1)(n−2)∫Sn−1D0f(α,θ)δ(n−2)(〈ω,θ〉)dS(θ),(12)∂(n−2)∂s(n−2)Rfα(ω,〈α,ω〉)=(−1)(n−2)∫Sn−1D1f(α,θ)δ(n−2)(〈ω,θ〉)dS(θ),
where α∈Γ,ω∈Sn−1,dS denotes the surface
measure on Sn−1,R is the n-dimensional Radon transform defined by(13)(Rf)(θ,s)=∫〈x,θ〉=sf(x)dx,and(14)fα(x)=〈f(x),|x−α|−1(x−α)〉ℝnis the projection of f onto (x−α)/|x−α|.In the following, we want to limit our
presentation ton=3.
Then formula (12) reads as
(15)∂∂sRfα(ω,〈α,ω〉)=(−1)∫S2D1f(α,θ)δ'(〈ω,θ〉)dS(θ)=∫S2∩{〈θ,ω〉=0}〈∇yD1f(α,y=θ),ω〉dS(θ).In the scalar case a solver forD0 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 Rg(ω,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 D0 since then f(x) is independent of α but unfortunately that does not help in case
of D1 (and analogous transforms for tensor fields as
well, see [13]), since
there the object function fα=〈f(x),(x−α)/|x−α|〉ℝ3 of R depends on and hence changes with α, see (15). Thus we seek an alternative way of solving D1f=g for vector fields f.


## 3. Approximation of Reconstruction Kernels  in Vector Field Tomography

As already said, an inversion formula forD1 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 D1. 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 D1. 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 ofapproximate 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 D0.
We describe this approach and then formulate an extension of it to D1.Letf∈L2(Ω3) be a scalar function. The approximate inverse
computes a smoothed version fγ of f by convolving f with a mollifiereγ∈𝒞∞(ℝ3). A mollifier eγ is a smooth function with small essential
support having the property that
(16)fγ(x):=(f∗eγ)(x)→f(x)asγ→0.Here, ∗ denotes the convolution(17)(f∗h)(x)=∫ℝ3f(y−x)h(y)dy.
Such a function is given by the Gaussian kernel
(18)eγ(x)=γ−3(2π)3/2exp(−|x|22γ2).
Provided that we can solve the equation
(19)D0∗[vγ(x)]=eγ(x−⋅),
then we can reconstruct fγ from the measured cone beam data D0f by
(20)fγ(x)=〈D0f,vγ(x)〉L2(Γ×S2)=∫Γ∫S2(D0f)(α,ω)vγ(x;α,ω)dS(ω)dα,
where vγ(x)=vγ(x;α,ω)∈L2(Γ×S2) for x∈Ωn is called a reconstruction kernel.
Hence the method of approximate inverse consists of evaluating inner products
of the given data D0f 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 toD1 and hence to VFT, we construct mollifier
fields Eγj∈L2(Ω3,ℝ3) defining
(21)Eγj(x):=eγ(x)⋅ej,j∈{1,2,3},
where e1=(1,0,0)⊤,e2=(0,1,0)⊤ and e3=(0,0,1)⊤. Using again the Gaussian (18) as mollifier eγ we obtain
(22)(fγ)j(x):=(f∗Eγj)(x)→fj(x)asγ→0
for f∈L2(Ω3,ℝ3). Unfortunately, by now the exact reconstruction kernels Vγj(x), that is, the solutions of D1∗[Vγj(x)]=Eγj(x−⋅) are still unknown. But the special structure
of the mollifier fields Eγj allow for a computation of reconstruction kernels for
(23)Pf(α,ω)=∫0∞f(α+tω)dt
with the help of kernels for D0.Theorem 3.
Letvγ be the reconstruction kernel associated to eγ with respect to D0, that is,
(24)D0∗[vγ(x)]=eγ(x−⋅).
Defining Vγj(x;α,ω):=vγ(x;α,ω)⋅ej∈L2(Γ×S2,ℝ3) yields
(25)P∗[Vγj(x)]=Eγj(x−⋅),
that means Vγj is a reconstruction kernel associated to Eγj with respect to P. The adjoint P∗ of P is given as
(26)P∗g(x)=∫Γ|x−α|−2g(α,x−α|x−α|)dα
for g∈L2(Γ×S2,ℝ3).Proof.
The adjointP∗ is computed as
(27)〈Pf,g〉L2(Γ×S2,ℝ3)=∫Γ∫S2∫0∞〈f(α+tω),g(α,ω)〉dtdS(ω)dα=∫Γ∫R3〈f(x),g(α,x−α|x−α|)〉|x−α|−2dxdα=〈f,P∗g〉L2(Ω,ℝ3),
where we applied Fubini’s
theorem as well as the substitution x=α+tω. A short calculation further shows that(28)P∗[Vγj(x)](y)=∫Γ|y−α|−2Vγj(x;α,y−α|y−α|)dα=∫Γ|y−α|−2vγ(x;α,y−α|y−α|)dα⋅ej=D0∗[vγ(x)](y)⋅ej=eγ(x−y)⋅ej=Eγj(x−y).The dataPf are not known and cannot be computed from D1f. But, observing that
(29)D1f(α,ω)=〈Pf(α,ω),ω〉ℝ3
and since Pf(α,ω)∈ℝ3 we obtain
(30)Pf(α,ω)=D1f(α,ω)ω+λ1(α,ω1⊥)ω1⊥+λ2(α,ω2⊥)ω2⊥,
where ω1⊥,ω2⊥∈S2 are such that {ω,ω1⊥,ω2⊥} is an orthonormal basis of ℝ3 and λ1,λ2 are appropriate coefficients. Thus
approximating
(31)Pf(α,ω)≈D1f(α,ω)ω
we neglect the parts orthogonal to ω and can apply the method of approximate
inverse using the reconstruction kernels Vγj for P.
This procedure results in Algorithm 1.Algorithm 1:For cone beam VFT.
Given : Measured data D1f(α,ω) for α∈Γ,ω∈S2Output    : Approximationfγ to fCompute:(1) g(α,ω)=D1f(α,ω)ω(2) (fγ)j(x)=〈g,Vγj(x)〉L2(Γ×S2,ℝ3)=∫Γ∫S2〈g(α,ω),Vγj(x;α,ω)〉dS(ω)dα=∫Γ∫S2D1f(α,ω)ωjvγ(x;α,ω)dS(ω)dαfor j∈{1,2,3}It is worth to mention that the mathematical model has
not been changed. Results for the scalar cone beam transformD0f(α,ω) are transferred to its 3D equivalent Pf(α,ω) which is an approximation to D1f(α,ω) by using Pf(α,ω)≈D1f(α,ω)ω.Figures3 and 4 display first results of the above
algorithm when applied to exact, simulated data for the solenoidal vector
fields f(x)=(−x2,x1,0)⊤ and f(x)=(1−x22−x32,0,0)⊤. As already said, we made use of the measurement setup as shown in Figure 1. The
scanning curve was Γ=rS2∩{x3=0},r=3, that is a circle of radius r=3 in the plane {x3=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 {x3=0}. The mollifier eγ defining the fields Eγj 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.(a) Original vector fieldf(x)=(−x2,x1,0)⊤. (b) Reconstruction with the described
algorithm for γ=0.00692 using exact, simulated data D1f.
(a)(b)(a) Original vector fieldf(x)=(1−x22−x32,0,0)⊤. (b) Reconstruction with the described
algorithm for γ=0.007 using exact, simulated data D1f.
(a)(b)Figure 5
Reconstruction kernelvγ(0;α,ω) for γ=0.007 associated to the scalar cone beam transform D0 according to [12].Figure4 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 x3-direction. Figures 6 and 7 show the
reconstructions for the planes {x3=0.5},{x3=0.75},{x3=0.9} and {x3=0.95} of the two afore-mentioned vector fields. Figure 7 also illustrates that the intensity
of the vector field f(x)=(1−x22−x32,0,0)⊤ decreases as should be expected since x32 is subtracted in the first component and that
even the directional error at the boundary is reduced with increasing x3.Reconstruction of the vector fieldf(x)=(−x2,x1,0)⊤ with the described algorithm for γ=0.00692 using exact, simulated data D1f. (a),(b) The planes {x3=0.5} and {x3=0.75}. (c),(d) The planes {x3=0.9} and {x3=0.95}.
(a)(b)(c)(d)Reconstruction of the vector fieldf(x)=(1−x22−x32,0,0)⊤ with the described algorithm for γ=0.007 using exact, simulated data D1f. (a),(b) The planes {x3=0.5} and {x3=0.75}. (c),(d) The planes {x3=0.9} and {x3=0.95}.
(a)(b)(c)(d)The very simple vector fieldf(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 {x3=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 x3-direction is shown in Figure 9. As in the reconstruction of the field f(x)=(1−x22−x32,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 {x3=0}. But the images also show that the intensity of the vector field is slightly
decreasing with increasing x3 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) Original vector fieldf(x)=(1,0,0)⊤. (b) Reconstruction with the described
algorithm for γ=0.0075 using exact, simulated data D1f.
(a)(b)Reconstruction of the vector fieldf(x)=(1,0,0)⊤ with the described algorithm for γ=0.0075 using exact, simulated data D1f. (a),(b) The planes {x3=0.5} and {x3=0.75}. (c),(d) The planes {x3=0.9} and {x3=0.95}.
(a)(b)(c)(d)Since the described algorithm allows the
reconstruction of any arbitrary plane in thex3-direction, a vertical cross-section of a
vector field can easily be computed. Figure 10 shows a reconstruction of the
plane {x1=0} for the circular vector field f(x)=(−x2,x1,0)⊤ for two different viewing angles. Figures 11
and 12 display the analogous results for the vector fields f(x)=(1−x22−x32,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.Reconstruction of the vector fieldf(x)=(−x2,x1,0)⊤ with the described algorithm for γ=0.00692 in the plane {x1=0}.
(a)(b)Reconstruction of the vector fieldf(x)=(1−x22−x32,0,0)⊤ with the described algorithm for γ=0.007 in the plane {x1=0}.
(a)(b)Reconstruction of the vector fieldf(x)=(1,0,0)⊤ with the described algorithm for γ=0.0075 in the plane {x1=0}.
(a)(b)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 {x3=0} we have to face another difficult problem.
Considering the plane {x3=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 ω∈{x3=0} as
(32)D1f(α,ω)=∫0∞〈f(α+tω),ω〉dt=∫0∞(f1(α+tω)ω1+f2(α+tω)ω2+f3(α+tω)ω3)dt=ω3=0∫0∞(f1(α+tω)ω1+f2(α+tω)ω2)dt,
which means that the vertical
component f3(x) of any vector field cannot be reconstructed in
that particular case. But even for planes where x3≠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 x1-x3-plane hitting the object Ω in just one point. From Figure 13 we easily see that
(33)sin(β)=Rr,tan(β)=ur⟺u=rtan(β)=rtan(arcsin(Rr)).
Using the equivalence
(34)arcsin(x)=arctan(x1−x2)
we obtain
(35)u=rtan(arctan(R/r1−(R/r)2))=R(r2−R2)/r2=rRr2−R2.
For the direction vector ω we can then deduce
(36)ω=1∥θ∥θwithθ=(−r0u),∥θ∥=r2+u2=r2+(rR)2r2−R2=r1+R2r2−R2=rr2−R2+R2r2−R2=r2r2−R2,ur2+u2=(rR/r2−R2)(r2/r2−R2)=Rr,−rr2+u2=−r(r2/r2−R2)=−(r2−R2)r⇒ω=(−(r2−R2)r0Rr).Figure 13
Clarifying the calculation of the maximal value forω3.This simple geometric calculation shows that the angleβ between the plane {x3=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
(37)−Rr≤ω3≤Rr⟺−13≤ω3≤13.
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 x3-direction as mentioned above.The consequences of the problem can easily be
illustrated by applying our algorithm to the two fully three-dimensional vector
fieldsf(x)=(1,1,1)⊤ and f(x)=(1,1,−1)⊤ whose x3-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 {x3=0}. Only for larger (or smaller) values of x3, that means planes above (or below) the plane of the scanning curve {x3=0}, the reconstructions become distinguishable. This is shown in Figure 15 where
the results for {x3=0.75} and {x3=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 x3-direction.(a),(b) Reconstruction of the vector fieldf(x)=(1,1,1)⊤ with the described algorithm for  γ=0.007 in the plane {x3=0}. (c),(d) Reconstruction of the vector field  f(x)=(1,1,−1)⊤ with the described algorithm for  γ=0.007 in the plane {x3=0}.
(a)(b)(c)(d)Top: reconstruction of the vector fieldf(x)=(1,1,1)⊤ with the described algorithm
for γ=0.007 in the planes (a) {x3=0.75}  and (b){x3=0.95}. Bottom: reconstruction of the vector field f(x)=(1,1,−1)⊤  with the described algorithm for γ=0.007
in the planes (c) {x3=0.75} and (d) {x3=0.95}.
(a)(b)(c)(d)Despite all these drawbacks it is nevertheless
possible to reconstruct a whole three-dimensional vector field. Summarizing the
results so far, thex1- and x2-component of any vector field f(x)=(f1,f2,f3)⊤, that is, f1 and f2, can be reconstructed for any plane in the x3-direction. Changing our measurement setup by
adding to the current scanning curve Γ=rS2∩{x3=0},r=3, a second, orthogonal circle of the same radius in the plane {x1=0} enables us to reconstruct the x2- and x3-component of any vector field, that is f2 and f3, 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 f1 and f2 for any plane in the x3-direction with the algorithm presented at the
beginning of the paper. The second, orthogonal circle in the plane {x1=0} enables us to analogously reconstruct f2 and f3 for any plane in the x1-direction. Taking a vertical cross-section at x3=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 f2-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 {x3=0} as usual, and one in the plane {x1=0}.(a) Reconstruction of the vector fieldf(x)=(1,1,1)⊤ with the described algorithm for γ=0.007 with circles in the planes {x3=0} and {x1=0}. (b) Reconstruction of the vector field f(x)=(1,1,−1)⊤ with the described algorithm for γ=0.007 with circles in the planes {x3=0} and {x1=0}.
(a)(b)It should be added that it is also possible to choose
the second orthogonal circle to lie in the plane{x2=0} instead of {x1=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) Reconstruction of the
vector fieldf(x)=(1,1,1)⊤ with the described algorithm for γ=0.007 with circles in the planes {x3=0} and {x2=0}. (b) Reconstruction of the vector field f(x)=(1,1,−1)⊤ with the described algorithm for γ=0.007 with circles in the planes {x3=0} and {x2=0}.
(a)(b)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 Figure19 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 {xi=0},i=1,2,3.(a) Reconstruction of the vector fieldf(x)=(1,1,1)⊤ with the described algorithm for γ=0.007 using three orthogonal circles in the planes {xi=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 {xi=0},i=1,2,3.
(a)(b)To verify our assumptions we reconstructed once again
the planar circular vector fieldf(x)=(−x2,x1,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 {x3=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) Original vector fieldf(x)=(−x2,x1,0)⊤. (b) Reconstruction with the described
algorithm for γ=0.00692 using one single circle in the plane {x3=0}. (c) Reconstruction with the described algorithm for γ=0.00692 using circles in the planes {x3=0} and {x1=0}. (d) Reconstruction with the described
algorithm for γ=0.00692 using circles in the planes {x3=0} and {x2=0}. (e) Reconstruction with the described
algorithm for γ=0.00692 using three orthogonal circles in the planes {xi=0},i=1,2,3.
(a)(b)(c)(d)(e)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 Dm. Denisjuk also proved in [30] that a full reconstruction of solenoidal vector
fields is possible if the vectorial cone beam data D1f(α,ω) are available for all α on a trajectory Γ satisfying a generalized Tuy condition and all
directions ω∈S2.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.



---

*Source: 174283-2009-02-03.xml*