Nonreciprocal broken ray transforms with applications to fluorescence imaging

Abstract Broken ray transforms (BRTs) are typically considered to be reciprocal, meaning that the transform is independent of the direction in which a photon travels along a given broken ray. However, if the photon can change its energy (or be absorbed and re-radiated at a different frequency) at the vertex of the ray, then reciprocity is lost. In optics, non-reciprocal BRTs are applicable to imaging problems with fluorescent contrast agents. In the case of x-ray imaging, problems with single Compton scattering also give rise to non-reciprocal BRTs. In this paper, we focus on tomographic optical fluorescence imaging and show that, by reversing the path of a photon and using the non-reciprocity of the data function, we can reconstruct simultaneously and independently all optical properties of the medium (the intrinsic attenuation coefficients at the excitation and the fluorescence frequency and the concentration of the contrast agent). Our results are also applicable to inverting BRTs that arise due to single Compton scattering.


Introduction
The broken-ray transform (BRT), also referred to as the V-line transform, has attracted significant recent attention [1][2][3][4][5][6][7][8][9]. This transform is a generalization of the classical Radon transform to the physical setting in which a photon can change its propagation direction due to scattering. A related recent development is inversion of the conical Radon transform, which arises in tomographic applications of the Compton camera [10][11][12][13][14][15][16][17][18], although in the latter case the vertex of the V-line is located at the detection surface rather than within the medium.
The BRT is defined in terms of integrals of a function (the attenuation coefficient) along two rays with a common vertex. The ray directions are fixed and determined by the directions of the source and detector. However, the position of the vertex can be scanned over the plane that contains both rays. The corresponding data provides sufficient information to reconstruct the attenuation coefficient in the plane thus defined, assuming that the scattering coefficient is spatially homogeneous. An interesting feature of the BRT is that its inversion does not require data from multiple projections; a single scan, analogous to one plane-parallel projection of traditional x-ray tomography, suffices for reconstruction [19]. In the single-projection geometry, inversion of the BRT is a mildly ill-posed problem [20].
If the scattering coefficient of the medium is not spatially homogeneous, more elaborate methods can be employed to reconstruct the scattering and the attenuation coefficients separately. This is possible in a measurement geometry with several (more than two) rays having a common vertex and lying in the same plane [21]. A direct and stable inversion formula for a transform involving four such rays has been derived in [22]. The formula is local, that is, it involves only first-order derivatives with respect to the vertex position, and therefore it does not require a 'complete' data set. This property is very useful and, in principle, is unattainable in transforms involving only integrals along lines. A generalization of the local inversion form ula to the case of an arbitrary number of rays, which defines the star transform, was derived in [23]. In particular, it was shown that a local formula can be obtained with only three rays. This imaging geometry is considered below. We note that a single 'star' consists of several distinct broken rays, some of which can share one edge. In a given measurement (utilizing a fixed source-detector pair), a photon always travels along one of the broken rays. We can utilize several detector arrays to perform measurements corresponding to different broken rays belonging to the same 'star' in one scan.
In previously considered applications, the BRT is reciprocal, which means that the transform is independent of the direction in which the photon travels along a given broken ray. Consequently, the associated measurements are independent of the interchange of the source and detector. However, if the photon is scattered or absorbed and re-emitted at a different frequency, the reciprocity of the BRT is lost. In this paper, we consider a non-reciprocal BRT applicable to the problem of fluorescence imaging in a weakly-scattering or non-scattering medium. In this case, the photon travels along a line until it is absorbed by a fluorophore molecule and then is re-emitted in a different direction and at a different (generally lower) frequency. We show that the non-reciprocity of the BRT allows us to access additional information about the medium by interchanging the source and detector in each pair. In particular, this approach allows us to reconstruct the attenuation coefficient of the medium at both the excitation and the fluorescence frequencies independently, as well as the concentration of the fluorophore. The method described in this paper requires doubling the number of physical measurements, but it does not involve any assumptions about the spectral dependence of the attenuation coefficient, which can generally differ at different points in the medium.
The proposed method is also applicable to x-ray imaging where the photon can change its direction due to Compton scattering. As is well known, the latter is accompanied by a reduction of the photon energy. Since the attenuation coefficient can depend on energy, this poses an additional challenge. This problem was solved in [24] under the assumption that the spectral dependence of the attenuation coefficient is linear. Here we show that accounting for the non-reciprocity of the measurements allows one to avoid making assumptions about the spectral dependence of the attenuation. However, in this case, we are restricted to using only imaging geometries in which the scattering angles (and hence the Compton energy shifts) are the same. The three-ray star geometry considered below satisfies this requirement.
The remainder of this paper is organized as follows. In section 2, we consider the coupled transport equations and show that, in a non-scattering or a weakly-scattering medium, the specific intensity of light at the fluorescence frequency due to a collimated incident beam at the excitation frequency is mathematically related to a BRT of the medium. In section 3, we define the data function to be used in the reconstruction of the functions of interest. In section 4, we consider a more general arrangement of sources and detectors in which several broken rays with a common vertex form a 'star'. In this section we show how the BRT non-reciprocity can be used to formulate the inverse problem. In section 5, we show how all three functions of interest (the background attenuation coefficients at the excitation and fluorescence frequencies and the concentration of fluorophore) can be reconstructed stably and separately by accounting for the BRT non-reciprocity. Explicit image reconstruction formulas are derived. In this section, we assume that the 'complete data' are available. In section 6, we consider physical situations in which the complete data are not available and sketch some approaches to performing partial or complete reconstructions in this less favorable case. Finally, section 7 contains a summary of the obtained results.

Coupled transport equations
We begin by considering the physical problem of fluorescence imaging in a weakly-scattering or non-scattering medium, taken to be a three-dimensional bounded domain Ω. In the absence of scattering, the transport equations describing the propagation of light at the excitation and fluorescent frequencies (distinguished by the subscripts e and f, respectively) are of the form Here I e (r,ŝ), I f (r,ŝ) are the specific intensities of light at the position r in the direction ŝ, and µ e , µ f are the attenuation coefficients of the medium at the frequencies indicated by the subscripts, and S is the fluorescent source. Equations (1a) and (1b) are also supplemented by the half-range boundary conditions I e (r,ŝ) = I inc (r,ŝ) and where n is the outward unit normal to ∂Ω and I inc is the incident specific intensity at the excitation frequency. Note that no light at the fluorescent frequency enters the medium except due to the source S, which is defined below. For a medium containing fluorescent molecules of the number density n(r), we can write where µ (0) e (r) and µ (0) f (r) are the intrinsic (background) attenuation coefficients of the medium and σ e , σ f are the absorption cross sections of the fluorescent molecules at the frequencies indicated by the subscripts. The cross sections are assumed to be known from spectroscopic measurements and it is expected that σ f < σ e . The source function at the fluorescence frequency is given by where η is the quantum efficiency of the fluorophores and the energy density u e (r) = I e (r,ŝ)d 2 s. If the incident light is injected into the medium as a narrow collimated beam, then it can be seen that [19] Here r 1 and s 1 are the location and direction of the source, W is the source power, and b a f d is the integral of f along the line connecting the two points a and b.
Let the specific intensity at the fluorescence frequency be measured at the point r 2 ∈ ∂Ω and in the direction ŝ 2 . Then it follows from (1b) that where we have accounted for the isotropy of the source function S. An integral similar to (6) has been computed in the appendix of [19]. Here we adduce the final result, Various geometrical quantities and objects appearing in (7) are illustrated in figure 1. In particular, sin(θ1+θ2)ŝ 2 is the position of the vertex, that is, the point where the two rays shown in figure 1 by the dashed blue lines intersect, Θ(·) is the unit step function, and we have assumed that the vectors r 21 and ŝ 1 lie in the XZ plane of the laboratory frame. The delta-function δ(ϕ 2 − π) ensures that ŝ 2 lies in the same plane. Therefore, equation (7) yields a nonzero result only if the three vectors r 21 , ŝ 1 and ŝ 2 lie in the same plane, which defines the slice of the medium in which the two-dimensional reconstruction of µ (0) e (r), µ (0) f (r) and n(r) is performed. Finally, the function Θ(π − θ 1 + θ 2 ) is equal to unity if the vertex exists and to zero otherwise. We note that, when experimental data are used, accounting for the geometrical factor 1/r 21 sin θ 1 sin θ 2 in (7) is important for correct formulation of the inverse problem.

Data function
Our goal is to reconstruct the three functions µ f (r) and n(r) separately from collimated boundary measurements of the fluorescence intensity. In an experiment, the latter can be registered with the use of a spectral filter that excludes all radiation at the excitation frequency. To proceed, it is convenient to introduce the dimensionless concentration of fluorophores according to ñ(r) = ησ 3/2 e n(r). We note that ñ ∼ 1 corresponds to a very high concentration and, in practice, we expect that ñ 1 will hold. Let us assume that the condition θ 1 + θ 2 < π holds so that the vertex exists. We then have Both sides of this equation are dimensionless. However, the expression is still singular and not amenable to direct interpretation as a measurable signal. To alleviate this problem, we note that the fluorescence intensity at the point of observation r 2 depends on the spherical coordinates of the unit vector ŝ 2 , θ 2 and φ 2 . The dependence on θ 2 is slow (1/ sin θ 2 ) while the dependence on ϕ 2 is fast and expressed mathematically by the delta-function δ(ϕ 2 − π). On the other hand, all physical detectors measure the specific intensity in some finite solid angle. We can assume therefore that the measured quantity is π+δ π−δ I f dϕ 2 , where δ is a small angle. We can now define the data function φ 12 (R) as follows: Let ŝ 1 and ŝ 2 be fixed and such that the vertex exists. Then there is a one-to-one correspondence between the position of the vertex R and the pair of variables (r 1 , r 2 ). In what follows, we will view R as the independent variable. Applying the definition (9) to (8), we obtain the equation where Here we have assumed that µ e (r) and µ f (r) are supported in Ω so that the ray integrals can be extended to infinity. Note that the data function can be measured only for such positions of the vertex R that ñ(R) = 0. If ñ(R) = 0, the measured fluorescence intensity is zero, at least, within the approximation used here, and the logarithm in (9) is undefined.

Inverse problem
Above, we have considered one broken ray defined by the locations of the source and detector (r 1 and r 2 ) and the location of the vertex r (from now on, we denote the vertex position by the small letter r). We now allow a more general arrangement of source-detector pairs, assuming that all broken rays corresponding to different source-detector pairs have the same vertex and that this vertex can be scanned with a reasonable spatial resolution over the sub-region of Ω in which ñ(r) > 0. Some of the broken rays in a given 'star' can have a common edge. The imaging geometry is sketched in figure 2, where û k are unit vectors pointing from the vertex to a source or detector. We note that û k = ±ŝ k , where the plus sign must be chosen if the ray points towards a detector and the minus sign is chosen if the ray points towards a source.
The generalization of (10) to the above case is the N-ray star transform equation Here Notably, the data function φ kl (r) is not symmetric with respect to the interchange of the indexes k and l. Physically, this means that the measurements do not obey reciprocity under the interchange of the source and the detector. The lack of reciprocity follows from the spectral dependence of the attenuation, that is, from µ e (r) = µ f (r). This is different from the star transform that arises in single-scattering tomography [23]. In the latter case, the data function is symmetric. In the case considered here, the lack of symmetry of the data function allows one to derive a local reconstruction algorithm, similar to those described in [22,23], which yields all three functions µ e (r), µ f (r) and ñ(r). Let us consider the symmetric and anti-symmetric linear combinations of the data, It should be kept in mind that measuring these linear combinations requires reversing the path of the photon, that is, physically interchanging the source and detector for each broken ray or, alternatively, performing an additional scan with interchanged source and detector arrays. From (12) and (14), we find the following set of linear equations: where

Inversion with complete data
As mentioned in section 3, the data functions φ kl (r) or φ (±) kl (r) are defined only for positions of the vertex r such that ñ(r) > 0. Otherwise, the data function can not be measured. In practice, the condition that the data is measurable is stronger and reads ñ(r) ñ min (r) > 0. Here ñ min (r) is experimentally determined.
The above considerations are illustrated in figure 3. Let Σ ⊂ Ω be the region in Ω in which the fluorescence signal can be measured. Σ can be defined with or without reference to ñ min (r); in an experiment, Σ is the set of all vertex positions for which the fluorescence signal is well above the noise level. In figure 3, Σ is shown as a simple connected region but one can consider more complex geometries. If Σ = Ω, the fluorescence signal can be measured in the whole domain of interest. We say that in this case we have access to complete data. This situation occurs when a significant concentration of the fluorescent molecules are present everywhere in Ω. In this section, we discuss inversion of (16) with complete data. Then in section 6 we sketch some approaches to reconstruction when complete data are not available. We begin by noting that (16a) is equivalent to the star transform [23] of a medium with the attenuation coefficient μ(r). The dimensionless concentration of fluorophores ñ(r) plays the role of the contrast of the scattering coefficient. We can therefore use the methods described in [23] to reconstruct simultaneously μ(r) and ξ(r); the dimensionless concentration is then trivially computed as ñ(r) = exp(−ξ(r)).
Consider the three-ray geometry shown in figure 2. We then define the vector coefficients a kl (k, l = 1, 2, 3) as shown in the matrix  Here the scalar coefficients σ k satisfy the condition Such coefficients can always be found since three vectors on a plane can not be linearly-independent. We then take the following linear combination of the symmetric functions φ (+) kl (r) (by symmetry, we mean here the property φ It is easy to see that the resulting function Φ (+) (r) satisfies We now use the property −∇ · u k I (+) k (r) =μ(r) to find that Thus, we have a local reconstruction formula for the spectrally-averaged attenuation coefficient μ = (µ e + µ f )/2. It is not surprising that this formula utilizes the symmetric data functions φ (+) kl (r) or a linear combination thereof Φ (+) (r), as defined in (20). Indeed, the symmetric data points are obtained by combining the measurements in which the photon travels through the medium in both possible directions (for a given broken ray). When the photon travels in one of these directions, it is attenuated at the rate µ e in the first segment of the broken ray and µ f in the second segment. When the direction of propagation is reversed, the rate of attenuation is µ f in the first segment and µ e in the second segment. By adding the two measurements together, we are effectively performing an average of the attenuation coefficient.
Since we have assumed here that complete data are available, we can reconstruct μ(r) everywhere in Ω. We can then use this result to compute the ray integrals I (+) k (r) according to (16a). Once this is done, we can obtain ξ(r) from any of the equations in (15a). This procedure yields a complete reconstruction of ñ(r). However, to obtain the spectrally-resolved intrinsic attenuation coefficients µ (0) e (r) and µ (0) f (r), we also need to know the spectral difference ∆(r) of the attenuation coefficients. We can use the anti-symmetric data functions φ (−) kl (r) to reconstruct ∆(r) by inverting (15b). The latter equation is similar to the star transform (15a) but has a different sign and also does not contain the term ξ(r). Therefore, inverting (15b) is slightly different from inverting (15a).
Recall that in order to derive a local inversion formula for (15a), we need to find a set of vector coefficients a kl that satisfy the following four conditions [23]: (i) a kl = a lk ; (ii) a kk = 0; (iii) l a kl = σ kûk ; and (iv) kl a kl = 0. For the conditions (i) and (iv) to be consistent, the scalar coefficients σ k must satisfy k σ kûk = 0. If the above four conditions hold, we have 1 2 kl a kl (x k + x l + y) = k σ kûk x k , where x k and y are arbitrary numbers. This property allows us to invert the star transform (15a) locally.
In order to invert (15b), we need to find a set of coefficients b kl such that It is easy to see that b kl must satisfy only three conditions, namely, (i) b kl = −b lk ; (ii) b kk = 0; and (iii) l b kl = σ kûk . We note that the condition (iv) kl b kl = 0 also holds but is a consequence of (i) rather than an independent condition. It also follows from (i) and (iii) that σ k must still satisfy k σ kûk = 0. For the three-ray geometry considered here, the following matrix of coefficients b kl satisfy all stated conditions: The coefficients σ k in this matrix are the same as in matrix (23) and determined from (19). We note that (23) does not contain û 2 explicitly but û 2 is not linearly independent of û 1 and û 3 . In fact, (23) is not a unique choice of coefficients; the conditions stated above for the elements b kl can be written equivalently as b 12 = b 21 + σ 2û2 and b 13 + b 23 = σ 3û3 . In deriving (23), we have made the choice b 13 = 0. However, we can also take b 23 = 0, which results in the matrix There are other possible choices in which all six off-diagonal elements of the matrix are non-zero and all three vectors û 1 , û 2 and û 3 are explicitly present. Noting this will restore the symmetry, which is seemingly lost in (23). Indeed, there is nothing special about the vector û 2 , which is not present explicitly in (23) or û 1 , which is not present in (24).
We now define a combination of the anti-symmetric data functions according to For the matrix of coefficients defined in (23), or for any other such matrix with any choice for the two indices, the right-hand side of (25) is of the form − with complex chemical and structural composition, this might not hold. Since different species of absorbers can have different optical spectra, the assumption µ (0) e (r)/µ (0) f (r) = κ may not be justified. In this case, additional spectral measurements could be employed to access the missing information. For instance, we have assumed that the light recorded by the detectors has passed through a spectral filter that eliminates light at the excitation frequency ω e . However, we can also employ spectrally-resolved detectors that can measure the light intensities at the excitation and fluorescence frequencies separately. Then, if the medium is weakly scattering, we can use intensity measurements at the excitation frequency to reconstruct µ e (r) everywhere in Ω. This does not require an additional scan, only a spectrally-resolved measurement of intensity.
Finally, we note that knowledge of µ e (r) everywhere in Ω is insufficient to reconstruct n(r). We also need to know µ f (r) in Ω. It appears that the only feasible approach to recover this function is to perform an additional scan using ω f as the excitation frequency. Fluorescence in this case is not excited and the signal at ω f is due to single scattering. If this additional scan is possible, then the ray integrals I (+) k (r) can be computed according to (16a) and we can further solve for ξ(r) just as in the case of complete data.

Summary
We have shown that the BRT is non-reciprocal if the photon can change its frequency at the vertex of the broken ray. The non-reciprocity provides an opportunity to gain additional information about the medium by interchanging the source and detector in each measurement. We note that the image reconstruction methods described in this paper are based on the original idea of Katsevich and Krylov [22], and are direct and spatially-local.
We have discussed applications of non-reciprocity to tomographic optical imaging of a medium containing a fluorescent contrast agent. In this case, the quantities of interest are the number density of the fluorophores and the intrinsic attenuation coefficients of the medium at the excitation and fluorescence frequencies. All three functions can be reconstructed simultaneously and independently by inverting the three-ray star transform of the medium (consisting of three broken rays with a common vertex and some common edges, as shown in figure 2), assuming that a complete data set is available. Physically, this condition means that there is a significant concentration of the fluorescent agent in the region of interest. If this condition is not satisfied, reconstructions can still be performed, as can a complete reconstruction with the use of additional spectrally-resolved measurements, as is discussed in section 6.
Note that the formalism presented above does not make any assumption about the angles of the broken rays (except if applied to Compton scattering of x-rays). As a result, both backscattering and transmission measurement geometries can be implemented. The back-scattering measurement geometry is important for in vivo imaging of mesocopic systems.
The method proposed in this paper is applicable to single scattering x-ray tomography in which the photon energy changes due to Compton scattering. In this case, complete data will almost always be available. In order to apply the reconstruction methods described in this paper to Compton scattering, it is important to use only broken rays with the same angle to guarantee that the Compton energy shifts are always the same. The symmetric three-ray star geometry satisfies this condition: all scattering angles in this case are equal to 2π/3. We note that the formalism presented here enables a new approach to dual-energy x-ray CT that could potentially overcome some of the limitations of this technique due to hardware, software, or dose constraints [25]. Indeed, for a Compton scattering angle of 2π/3, the energy separation between the primary (excitation) and scattered radiation is helpful for improved tissue differentiation based on simultaneous reconstruction of tissue properties at two energies. At the same time, this method can be realized with single-energy exposure and single-energy detection. Moreover, the spectral distribution of the detected scattered radiation can be much narrower than that of the detected lower-energy radiation conventionally used in dual-energy CT. This could be useful for ameliorating the artifacts associated with beam hardening.