Unique Tomographic Reconstruction of Vector Fields Using Boundary Data

The problem of reconstructing a vector field v(r) from its line integrals (through some domain D) is generally undetermined since v(r) is defined by two component functions. When v(r) is decomposed into its irrotational and solenoidal components, it is shown that the solenoidal part is uniquely determined by the line integrals of v(r). This is demonstrated in a particularly simple manner in the Fourier domain using a vector analog of the well-known projection slice theorem. In addition, under the constraint that v (r) is divergenceless in D, a formula for the scalar potential phi(r) is given in terms of the normal component of v(r) on the boundary D. An important application of vector tomography, i.e., a fluid velocity field from reciprocal acoustic travel time measurements or Doppler backscattering measurements, is considered.


I. INTRODUCTION HE conventional tomography problem may be de-
T scribed as the task of reconstructing an unknown scalar function f ( r ) from its projections, or equivalently, from a complete set of line integrals through some bounded domain off.The solution to this problem has resulted in a wealth of applications in many disciplines, such as x-ray and NMR tomography in biomedicine, and acoustic and seismic tomography in oceanography and geophysics.The related problem of reconstructing a vector field from its line integrals has received far less attention, not for lack of potential applications, but because the vector reconstruction problem has generally been regarded as underdetermined.This seems clear from the fact that a scalar function is determined uniquely by its Radon transform (its projections), whereas a 2-D vector field requires two component functions to be determined.
A particularly important application of the vector tomography problem is the mapping of fluid flow from reciprocal acoustic travel time measurements through some bounded region of the flow [11][12][13][14].Here the vector field the latter problem, acoustic signals are propagated in opposite directions and the difference in travel time is recorded.This travel time difference is sensitive to fluid motion, but is relatively insensitive to stationary velocity inhomogeneities in the medium.One can also conceive of the optical analog of this problem, in which changes in the optical path length of a laser beam directed through a region of flow are measured interferometrically.A third example of a vector tomography problem arises from the measurement of Doppler shifted waves scattered from particles carried along with the moving fluid.Formulations of the fluid-flow reconstruction problem based on travel time and Doppler measurements are developed in more detail in the next section.
Each of the above problems can be reduced to solving an equation, of the type (1) below, as shown in the following section.However, in previous attempts at mapping fluid flow based on travel time tomography, conventional (scalar) tomographic theory has invariably been applied [ 11- [3]; this leads to an underdetermined problem and a nonunique solution unless the flow field is known in advance to be of a particularly simple form.
To define the problem, suppose the domain of the unknown vector field is D and its boundary is aD.'We show that a 2-D vector field v(r) can be uniquely recovered given two types of information: a) the line integrals of u(r) through D , and b) the normal component of v(r) on the boundary aD, with the constraint that v(r) is divergenceless in D. In the fluid-flow reconstruction problem, the latter constraint will hold when no sources or sinks reside in D.* The key to the reconstruction problem is to decompose v into is nonrotational and solenoidal com- We then show that the solenoidal part V x Y is uniquely determined by the line integrals of v along all lines joining pairs of points on the boundary aD.This fact has been previously demonstrated, [4], [5] but here we use a Fourier domain derivation that is particularly simple.In this derivation, we obtain a vector analog of the conventional (scalar) projection slice theorem.The projection slice theorem has tra-

~~ ~~
to be reconstructed represents the fluid field' In 'do is not a physical boundary, but merely the locus of points between which the integration paths are defined.In the acoustic travel time problem, tJD is the closed curve on which the sources and receivers reside.'This is a consequence of the continuity equation which, in the absence of sources and sinks, reads V .( p v ) = 0 where p is density.In almost all cases of practical interests, p may be regarded as effectively constant (this is always true for liquids), thus implying V .v = 0. and a and b range over the boundary dD of D. In ( l ) , .i is a unit vector along the line joining the points a and b and dl is an element of path length along this line.For convenience, we assume the domain D is convex so that the integration paths intercept only the interior of D. Equation (1) states that T(a, b) is the path integral of the component of v along the line joining a and b.
The solution to (1) is facilitated by decomposing the field v(r) into its irrotational and solenoidal components (known as Helmholtz's theorem [6]): v(r) = V 9 + V x Y, where 9 ( r ) and Y (r) are scalar and vector potentials.
We show in the following that (1) is, in fact, insufficient to determine v uniquely; in particular, the tomographic data T(a, b) alone do not determine the irrotational component VCP.However, for arbitrary v , a complete set of line integrals does determine the vector potential Y (r) uniquely, and hence, the solenoidal component V X Y of v(r).We then show that, if v(r) is divergenceless (V v = 0) in the domain of reconstruction D , an explicit formula for the scalar potential 9 may be derived in terms of v on the boundary of D.
In a recent work by Braun and Hauck [5] the vector tomography problem is also considered, but the treatment here differs from their analysis in several basic respects.These authors demonstrate, as we do, that the solenoidal component of the vector field is determined from line integral data.They go on to show that the remaining component (corresponding to the irrotational component V 9 in this paper) can be recovered if a so-called "transverse measurement" can be performed that measures the pathintegral of the component of the vector field normal to the path.Unfortunately, this transverse component of a fluid flow field cannot be measured by travel time or Doppler techniques, and is thus unavailable in flow reconstruction problems using acoustic waves.3The present treatment shows that the "transverse measurement" is not needed if the field can be measured on the boundary and if the 'Braun and H a w k do note that the line integral of the transverse component can be measured using an optical Schlieren technique, but this is only practical in rather specialized problems having optical access.field is divergenceless in D ; in particular, the missing irrotational component can be recovered from boundary information alone (assuming a divergenceless field in D).
Another difference between the present treatment and that of Braun and Hauck is that these authors employ a Radon transform formulation, whereas the Fourier domain derivation presented here is noticeably simpler.

EXAMPLES OF VECTOR TOMOGRAPHY
We begin by showing that the problem of reconstructing a 2-D fluid flow field from acoustic travel time measurements reduces to solving an equation like (1).Let u(r) denote the fluid velocity field and let c(r) denote the possibly spatially varying sound speed.Suppose a sound pulse is propagated along a line between a source at rl and a receiver at r2, where rl and r2 are on the boundary dD of the region D. We assume that IuI << c everywhere in D and that the variations in the sound speed c(r) are sufficiently small and/or the path lengths are sufficiently short so that ray refraction over all paths may be neg l e ~t e d .~ The travel time, TI,, of a pulse along a path between rl and r2 can be expressed as the path integral where .i denotes the unit tangent vector along the acoustic path and dl is an element of path length.This equation neglects an error that is second order in IuI /c.Reference [ l ] contains the derivation of (2), though it should seem intuitively reasonable that (2) holds for small IuI / e since the denominator of the integrand of (2) expresses the local sound velocity in the moving medium as the sum of the sound speed, c(r), in the absence of flow and the component of the flow velocity along the acoustic path (i.e., Again, neglecting terms which are second order in u(r) ..i).
\ U \ / e , (2) can be expanded to yield Now, let TZ1 denote the travel time obtained by transmitting in the opposite direction (i.e., by interchanging source and receiver).In this case, .i changes sign, and from (3) Equation ( 5) is of the form (1) when we identify v(r) with -2u(r) /c(r)'.Assuming measurements are taken be-4Refraction of the paths can always be neglected for sufficiently small velocity inhomogeneities (or, in the case of flow, for small u / c , where u and c are, respectively, the flow and wave speeds) since, under the assumption of a straight path, the resulting error in the path integral due to refraction is second order in the inhomogeneity.For further discussion, see 171.tween all points rl and r2 on the boundary, it is worth pointing out that (4) could be used to recover the sound speed c(r) by means of conventional (scalar) tomography.Equation (4) also shows that, to first order in IuI / c , the sum of the reciprocal travel times is independent of fluid movement.
One can conceive of the optical analog of the previous problem in which the quantity measured is the change in optical path length of a collimated laser beam directed through the region of flow.Perturbations in the optical path length due to flow will manifest themselves as optical phase shifts, and such phase shifts can be measured interferometrically.
An approach based on the Doppler effect also reduces to a problem in vector tomography.This technique, which appears not to have been previously reported, again gives rise to a system of equations (l), as we show in the following.In the Doppler approach, both optical and acoustic versions can be imagined.To describe the technique, consider a collimated beam (acoustic or optical) directed through a region of moving fluid in which particulate matter is present.The purpose of the particulate matter is to provide a distribution of scattering centers which moves with the fluid and from which Doppler-shifted waves are scattered.We shall assume, for simplicity, that the density of scatterers is uniform.The scattered waves from particles intercepted by the collimated beam arrive back at the source and are detected there.Suppose the source is monochromatic with time dependence exp ( -i d ) and assume that this time dependence has been removed from the received signal by demodulation5 (or, in the optical case, by a suitable optical heterodyning operation).What remains is the relatively slowly varying Doppler-modulated signal.If the collimated beam points in the direction S, the Doppler shift suffered by a wave backscattered from a particle in the beam is 2wuL/c0, where uL = v S is the component of particle velocity along the beam direction, w is the frequency of the incident wave, and c0 is the wave speed in the stationary medium.The total Doppler signal then consists of the sum of contributions from waves scattered by the particles lying along the path of the collimated beam.The result can be expressed as the path integral

rc/L(t) = R iL exp [ i a ( r ) ] exp [i2wtv(r) * ? / c o l dl. (6)
In (6), L denotes the path along which the collimated beam is directed, v(r) is the spatially varying fluid velocity, R is a constant related to the average scattering cross section of the particles intercepted by the beam,6 and a ( r ) is a random phase factor.Although the beam has finite thick-'We assume that, in the demodulation operation, both the in-phase and quadrature components of the signal are obtained, so that the phase of the received signals is well defined, and hence, the exponential notation in (6)-(9) is meaningful.
?n the most general situation R could also vary with position if the particle density is not uniform, or if the beamwidth diverges noticeably.In this case, R would become part of the integrand of (6).
ness, (6) is written, for simplicity, as a 1-D path integral.This is a reasonable approximation if the beamwidth is assumed much smaller than the scale of the inhomogeneities within the medium to be reconstructed.It should be borne in mind that (6) contains an implicit integration over the beamwidth (not shown) that has no significant effect if the above assumption holds.The spatially varying phase factor a ( r ) in (6) is introduced to account for any phase shift suffered by the wave due to scattering from the particle at r, as well as the phase shift undergone due to propagation of the wave between source and scatterer and back to the detector.Finally, we assume that the scattering is incoherent, i.e., that the particles scatter with random and uncorrelated phase.As a consequence of the latter assumption, the phase shift a ( r ) is a random function of r.Now, differentiating (6) with respect to time brings down a factor proportional to v(r) S:

exp [i2wtv(r) ? / c O ] ( v ( r ) S) dl. (7)
Next, multiplying (7) by the complex conjugate of (6), (8) The assumption that the phase shifts due to scattering at different locations are uncorrelated implies that, on carrying out the ensemble average: (exp [ i ( a ( r ) -a ( r ' ) ) ] ) oc 6(rr ' ) where 6 ( . ) is the 3-D Dirac delta function.Using this in (€9, the double integral collapses to a single integral and the exponential disappears, leaving where K is a constant.Equation (9) again is of the form (1).
It is worth pointing out that flow imaging can also be accomplished by the more traditional approach of pulsed Doppler imaging, in which spatial resolution along the path of the beam is achieved by recording the time of arrival of a backscattered pulse.An inherent property of the pulsed Doppler approach, which limits the sensitivity of the flow measurement, is the ever present range Doppler ambiguity [ 8 ] .The tomographic approach to flow imaging, although requiring multiple intersecting paths, does not suffer from the ambiguity problem since here the ex-citation is assumed to be monochromatic.Thus theoretically, the tomographic method has the potential for achieving much higher velocity sensitivity than traditional pulsed Doppler imaging techniques.

IV. THE RECONSTRUCTION PROBLEM
Let us parameterize the line through the points a and b by its distance p from the origin and the unit vector A normal to the line, as illustrated in Fig. 1; this line is denoted by L( p , A), so (1) may be written as where .i is a unit vector pointing along the line L (and normal to A).It is convenient to set v(r) = 0 outside of the measurement domain D ; this is permitted since the line integrals are only assumed given through the interior of D. This has several purposes.It is now permissible to (mathematically) extend the integration path in (10) along the line L( p , ri) from -a to 00.In addition, we apply the Helmholtz decomposition v = V+ + V X Y to the mathematically truncated field (v zero outside D).This satisfies the requirement that v vanish at infinity which is necessary for the uniqueness of the Helmholtz decomposition.' It is useful to recall the well-known projection slice theorem at this stage.We then derive an analogous theorem for the vector problem.The projection slice theorem forms the theoretical basis of reconstruction algorithms in conventional (scalar) tomography.The theorem can be stated briefly as follows.Assumef(r) is the (scalar) function to be reconstructed.Then the values of the 2-D Fourier transform o f f ( r ) on lines through the origin in the 2-D Fourier domain are given by the 1-D Fourier transforms of the projections of f ( r ) (i.e., the 1-D Fourier transforms of T( p , A) computed with respect to p ) .To state this precisely, consider the line integral of the scalar functionf(r) along the path L( p , A): 3   which can be equivalently written where 6 ( e ) is the 1-D Dirac delta function.The argument of the delta function is zero for points that lie on the line 'The uniqueness of the Helmholtz's decomposition holds at any point where U is continuous and differentiable, and under the assumption that U vanishes at infinity [ 6 ] .Differentiability will hold for any physical flow field, and since D is bounded, ZJ vanishes at infinity, and thus uniqueness is guaranteed inside D. Uniqueness breaks down only for those points on the boundary 6D since the truncated field is discontinuous there.This is of no consequence.however.since uniqueness does hold inside D . the region where Y is to be reconstructed.Indeed.the real (i.e., physical) field (as opposed to the truncated field) is assumed to be accessible to measurement on the boundary, since this is where the sources and receivers reside.We also require the normal component of e on the boundary to compute the scalar potential by means of (21).L( p , A), as required.Now letf(k) denote the 2-D Fourier transform off(r), i.e.: f ( k ) = j j'f(r) exp (-ik r) d2r (12) where k is a vector in the 2-D Fourier domain.Then Fourier transforming T( p , A) in (1 l ) with respect to p gives T(k, A) = ss f ( r ) exp (-&A r) d 2 r and on comparing to (12) shows that (13) This is the traditional form of the projection slice theorem.Equation (13) states that, by Fourier transforming each projection T( p , A) with respect to p and letting the direction of views A range over 180", one can recover the 2-D Fourier transformf(k), and hence, obtainf(r) via an inverse 2-D Fourier transformation.
To obtain a vector field analogue of the projection slice theorem, we write (10) in a form similar to (1 1): Taking the 1-D Fourier transform of (14) with respect to p gives T(k, A) = .i * E(kA) (15)   where is the 2-D Fourier transform of v(r).We now employ the decomposition v(r) = V 9 + V x Y, and note that in two dimensions, the vector potential Y has a single compo-  16) into (15), noting that ri x i = .i and recalling .i * ri = 0, (15) reduces to Note, significantly, that 6 ( k ) drops out in the derivation of (17), so the Fourier transform of the vector potential ' P, and hence, the solenoidal component V X Y of v(r), can be recovered from the l-D Fourier transform of the projections F(k, ri) independently of the irrotational part VIP.Equation (17) is the vector version of the projection slice theorem, i.e., the vector analog of (13).
Our remaining task is to recover the scalar potential 9.
Under the constraint that V * v = 0 in D, it follows that 9 solves Laplace's equation ( V 2 9 = 0) in D, and hence, can be recovered from the boundary values of v on aD.
To obtain an explicit formula for 9 in terms of U on aD, we note that the potential functions 9 and Y can be expressed in terms of U as follows [6]: (2 1) can be used to reconstruct a vector field, we consider in the next section a simple example that can be treated analytically.
IV. AN ANALYTICAL EXAMPLE Consider the problem of a constant vector field representing steady flow in one direction.Though this problem is physically trivial, it illustrates nicely the application of ( 17) and ( 21).(In a realistic problem involving data, ( 17) and ( 21) would, of course, be evaluated by numerical means.)Suppose the field is v(r) = voa, and vo is a constant and a a unit vector.We consider a bounded domain D in which the line integrals of v through D are known and the boundary values of v on aD are also known.We assume nothing is known about U in the interior of D, and our aim is to reconstruct v in D from the above information.
Before proceeding it is worth re-emphasizing the need to restrict ourselves to a bounded region of reconstruction D. This is particularly evident here, since the Helmholtz theorem clearly fails for an unbounded constant vector field [6].However, since the exterior of D has no influence on the measurements, we can ignore the exterior and deJne v(r) = 0 there.For the (mathematically) truncated field ( U zero outside D), both the solenoidal and irrotational parts of v are well defined, and (17) and (2 1) can be used to reconstruct both components, as we now show.i o ;

T ( P , f i ) =
where A x 2 = .i, 2 being the unit vector normal to the plane of flow.T( p , 2 ) can be thought of as our tomo- performing an ensemble average (denoted by ( -)) over all scatterers, results in LL exp [i2wt(v(r) -v ( r ' ) ) .;/col -( ~( r ) * S) dl d l ' .

tFig. 1 .
Fig.1.The integration path L is parameterized by its distance p from the origin and by its direction i.
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 1, NO. 3, JULY 1992 nent; that is, Y = \ki, where 2 is the unit vector normal to the plane of reconstruction.Our objective is to recover both functions 9 (r) and \k (r), which together determine u(r).The 2-D Fourier transform of v(r) = V 9 + V X Y is aD Finally assuming that V in (20) vanishes, leaving v = 0 in D, the second integral 9 ( r ) = ($ g ( r I r ' ) v ( r ' ) ri ds (21) fi(k) = ik&(k) + i(k x z^)G(k) (16) where 6 (k) and $ (k) are, respectively, the 2-D Fourier transforms of the scalar functions 9 ( r ) and !P(r).Now substituting ( where U(r) = ss v ( r ' ) g ( r I r ' ) d 2 r r (19) and g ( r ( r ' ) = -In (Irr' 1 ) / 2 ~ is the Green's function for the 2-D Laplace's equation.That is, g solves D V 2 g ( r ( r ' ) = -6(rr ' ) where V 2 is the 2-D Laplacian and 6 (rr ' ) is the 2-D Dirac delta function.Inserting (19) in the first equation of (18), gives 9 ( r ) = -v .U(r) = -s i v ( r ' ) .Vg(r1r') d 2 r ' the first step, Vg = -V ' g was used, the prime indicating differentiating with respect to r ' ; the last step * v and employing the divergence theorem.In (20), A denotes the outward unit vector normal to the curve aD.follows on using the identity v V ' g = V ' -( g v )g v which gives the scalar potential 9 explicitly in terms of the normal component of v on the boundary aD.In summary, the solenoidal part of an arbitrary 2-D vector field is determined uniquely by a complete set of line integrals of v through D. In particular, the 2-D Fourier transform of the solenoidal component, \k(kA), is given in terms of the l-D Fourier transforms of the projections F(k, A ) , by (17).If we impose the additional constraint that the field is divergenceless in D, then the irrotational part can be recovered, using (21), from the normal component of v on the boundary aD.As noted earlier, the restriction V v = 0 holds in the absence of sources and sinks in D. As an illustration of how (17) and For simplicity, let the boundary aD be a circle of radius R centered at the origin.Substituting v(r) = voa into (10) and computing the line integrals through D, gives 2vo(a .i>m,for I p~ < R for I pI > R (22)