The Broken Ray Transform: Additional Properties and New Inversion Formula

TL;DR

This paper explores the mathematical properties of the broken ray transform (BRT), introduces new inversion formulas leveraging Fourier analysis, and discusses practical implications for imaging modalities like optical and x-ray systems.

Contribution

It provides a new perspective on the inverse problem of BRT, compares existing inversion formulas, and derives computationally efficient inversion methods for arbitrary scatter angles.

Findings

New inversion formulas derived using Fourier transform

Numerical simulations demonstrate practical effectiveness

Clarification of data requirements for global reconstruction

Abstract

The significance of the broken ray transform (BRT) is due to its occurrence in a number of modalities spanning optical, x-ray, and nuclear imaging. When data are indexed by the scatter location, the BRT is both linear and shift invariant. Analyzing the BRT as a linear system provides a new perspective on the inverse problem. In this framework we contrast prior inversion formulas and identify numerical issues. This has practical benefits as well. We clarify the extent of data required for global reconstruction by decomposing the BRT as a linear combination of cone beam transforms. Additionally we leverage the two dimensional Fourier transform to derive new inversion formulas that are computationally efficient for arbitrary scatter angles. Results of numerical simulations are presented.

Equations

$$\exp \left(-{\int }_0^{\infty }\mu (x-t{\theta }_i)dt\right).$$

$$q(s,E)=2\frac{E}{hc}\sqrt{\frac{1-s}{2}}.$$

$$p(x,{\theta }_i,{\theta }_j,E)=f\left(x,q\left(-{\theta }_i\cdot {\theta }_j,E\right)\right)\exp \left(-{\int }_0^{\infty }\mu (x+t{\theta }_i)+\mu (x+t{\theta }_j)dt\right).$$

$$(B\mu )(x,\theta )\text{mathchar}58={\int }_0^{\infty }\mu (x+t\theta )dt.$$

$$(G\mu )(x,{\theta }_i,{\theta }_j)\text{mathchar}58=(B\mu )(x,{\theta }_i)+(B\mu )(x,{\theta }_j).$$

$$(R\mu )(v,\theta )\text{mathchar}58={\int }_{-\infty }^{\infty }\mu (v{\theta }^{\perp }+t\theta )dt.$$

$$\ln p(x,{\theta }_i,{\theta }_j,E)$$

$$=\ln f\left(x,q\left({\theta }_i\cdot {\theta }_j,E\right)\right)-(G\mu )(x,{\theta }_i,{\theta }_j).$$

$${\theta }_i\cdot {\theta }_k={\theta }_j\cdot {\theta }_k,$$

$$\ln p(x,{\theta }_i,{\theta }_k,E)-\ln p(x,{\theta }_j,{\theta }_k,E)=-(B\mu )(x,{\theta }_i)+(B\mu )(x,{\theta }_j).$$

$$\mu_{C}(x)=\left\{\begin{array}[]{ll}\mu(x),&\mbox{for }x\in C\\ 0,&\mbox{otherwise.}\end{array}\right.$$

$${\hat{\hat{b}}}_{\theta }(w)$$

$$={\hat{\hat{\mu }}}_C(w)\left[\frac{-1}{i2\pi w\cdot \theta }+\frac{1}{2}\delta (w\cdot \theta )\right].$$

$${\hat{\hat{g}}}_{i,j}(w)$$

$$={\mathcal{F}}^2\left\{(B{\mu }_C)(x,{\theta }_i)+(B{\mu }_C)(x,{\theta }_j)\right\}$$

$$={\hat{\hat{\mu }}}_C(w)\left[\frac{-w\cdot \left({\theta }_i+{\theta }_j\right)}{i2\pi \left(w\cdot {\theta }_i\right)\left(w\cdot {\theta }_j\right)}+\frac{1}{2}\delta \left(w\cdot {\theta }_i\right)+\frac{1}{2}\delta \left(w\cdot {\theta }_j\right)\right].$$

$${\hat{\hat{h}}}_{i,j}(w)\text{mathchar}58=\frac{-w\cdot \left({\theta }_i+{\theta }_j\right)}{i2\pi \left(w\cdot {\theta }_i\right)\left(w\cdot {\theta }_j\right)}.$$

$${\Theta }_{i,j}\text{mathchar}58=\left\{w\text{mathchar}58w\cdot \left({\theta }_i+{\theta }_j\right)=0,w\cdot w>0\right\}.$$

$${\hat{\hat{\mu }}}_C(w)={\hat{\hat{g}}}_{i,j}(w)\frac{-i2\pi \left(w\cdot {\theta }_i\right)\left(w\cdot {\theta }_j\right)}{w\cdot \left({\theta }_i+{\theta }_j\right)},\forall w\notin {\Theta }_{i,j}.$$

$${\mu }_C(x)$$

$$=\frac{-1}{\Vert {\theta }_i+{\theta }_j\Vert }\frac{d}{d{\theta }_i}\frac{d}{d{\theta }_j}{\int }_{-\infty }^0{g}_{i,j}\left(x+s\frac{\left({\theta }_i+{\theta }_j\right)}{\Vert {\theta }_i+{\theta }_j\Vert }\right)ds.$$

$$v_{\theta }^{-}$$

$$v_{\theta }^{+}$$

$$v_{\theta }$$

$$u_{\theta }^{-}(v)\text{mathchar}58=\min t\text{, s.t. }t\theta +v{\theta }^{\perp }\in C$$

$$u_{\theta }^{+}(v)\text{mathchar}58=\max t\text{, s.t. }t\theta +v{\theta }^{\perp }\in C.$$

$$f_{\theta }^{-}(v)\text{mathchar}58=u_{\theta }^{-}(v)\theta +v{\theta }^{\perp }$$

$$f_{\theta }^{+}(v)\text{mathchar}58=u_{\theta }^{+}(v)\theta +v{\theta }^{\perp }.$$

$$C_{\theta }^{-}$$

$$C_{\theta }^{+}$$

$$V_{\theta }^{-}$$

$$V_{\theta }^{+}$$

$$\overline{C\cup C_{\theta }^{-}}=C_{\theta }^{+}\cup V_{\theta }^{+}\cup V_{\theta }^{-}.$$

$$(B{\mu }_C)(x,\theta )={\int }_{-\infty }^0{\mu }_C(x+t\theta )+{\int }_0^{\infty }{\mu }_C(x+t\theta )dt$$

$$(B{\mu }_C)(x,\theta )$$

$$=(R{\mu }_C)(x\cdot {\theta }^{\perp },\theta ).$$

$$(R{\mu }_C)(v,\theta )$$

Full text

1 Introduction

The broken ray transform (BRT) appears in the forward model of a number of imaging modalities and measurement geometries. It was first considered in the context of optical scatter imaging [1] and later applied to x-ray scatter imaging [2]. The BRT occurs when the measured data are characterized by two rays sharing a common vertex. Time-of-flight positron emission tomography (TOF-PET) could also be considered in this framework (subject to the TOF ambiguity profile). The BRT has been considered for both translation-only measurement geometries [1, 3, 4, 2, 5, 6, 7, 8] and a rotate-shift measurement geometry [9, 7, 10].

To our knowledge, all applications associated with the single-scatter BRT represent joint reconstruction problems. Two spatially-varying images must be resolved. For example, it may be necessary to recover attenuation despite nonuniform scatter density [3, 2], or recover two attenuation images at distinct energy levels [11]. Variations in the joint reconstruction problem have motivated several novel contributions related to the BRT. These contributions are not strictly academic. In application, the forward models must be tailored to the joint reconstruction problem. This may limit available or useful data. The joint reconstruction problem determines available BRT inversion strategies.

Before contrasting prior work, we first define a notional measurement geometry. We will use this to establish some notation and define a joint reconstruction problem involving the BRT. We generalize the joint reconstruction problem to cover coherent scatter x-ray imaging. This modality has received renewed interest recently, however joint reconstruction of scatter density and attenuation has not yet been addressed [12, 13, 14]. In coherent scatter x-ray imaging, scatter density is highly sensitive to scatter angle. For this reason, we focus on BRT inversion using only two scatter angles.

As a simplification we assume a mono-chromatic x-ray pencil-beam incident upon some media of interest. At point ${x\in\mathbb{R}^{2}}$ the beam interacts with the media and scatters coherently. We use $\theta_{i}\in S^{1}$ to represent the direction of the source relative to the scatter location. The direction of the scattered photon is ${\theta_{j}\in S^{1}}$. We assume it is detected by a columated detector. Due to the combination of a pencil beam and columated detector we assume the scatter location $x$ is known precisely. This measurement geometry is depicted in Figure 1a.

The intensity measured at the detector largely depends on two media-specific images: attenuation and scatter density. The incident path is a straight line defined by $\theta_{i},x$. The loss in intensity along this path due to attenuation is governed by Beer’s law

[TABLE]

We use $\mu(x)\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{2}\rightarrow\mathbb{R}^{+}$ as the attenuation image representing both scatter and absorption. Intensity loss along the scatter path due to attenuation has a similar form and combines multiplicatively. For non-coherent scatter applications (e.g. fluorescence imaging) it may be necessary to distinguish the energy levels of the attenuation image before and after the scatter event. This has been investigated recently [11].

Even in homogeneous media, the intensity observed at the detector may vary with respect to scatter angle (e.g. $\theta_{i}\cdot\theta_{j}$) and energy level. For coherent scatter imaging, the scatter density does not depend on these terms independently, but rather through Bragg’s law [15]. This relationship is summarized by the so-called momentum transfer ${q(s,E)\mathrel{\mathop{\mathchar 58\relax}}(-1,1)\times\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}}$

[TABLE]

Here $h$ and $c$ are the Plank’s constant and the speed of light, respectively. This definition is unconventional as we have chosen to define it over the cosine of the scatter angle, $s$, rather than the scatter angle directly. Scatter intensity for inhomogeneous media varies both spatially and with respect to momentum transfer. We use $f(x,q)\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{2}\times\mathbb{R}^{+}\rightarrow\mathbb{R}^{+}$ to represent the scatter density image.

Combining the effects of attenuation and scatter density we arrive at the measurement function

[TABLE]

In this expression we have omitted a number of terms necessary for accurate models of measured data. However, we assume the remaining terms are known multiplicative factors. Measured data can then be scaled to achieve this generalized form.

To simplify the notation we will make use of three common transforms. Borrowing the notation of Natterer [16], we define the cone beam transform (CBT) $B$ of $\mu$

[TABLE]

This transform appears in (3). In particular, the generalized model includes the linear combination of two cone beam transforms sharing a common vertex. This is commonly referred to as the broken ray transform

[TABLE]

Denoting the left side of (5) by data $g_{i,j}$, Figure 1b represents (5) as a linear system with component CBT operators (4).

In subsequent sections we will also make use of the 2D Radon transform

[TABLE]

Here $v\in\mathbb{R}^{1}$, and $\theta\in S^{1}$ represent the shift and rotate coordinates of the transform. We assume $\theta^{\perp}$ is uniquely defined by rotating $\theta$ counter-clockwise by $\pi/2$.

Using these transforms we can express the log of the measured data

[TABLE]

The BRT is not directly available in (7b). However, the term $f$ can be canceled with differential measurements [3] even for inhomogenous media.

Given three scatter angles $\theta_{i},\theta_{j},\theta_{k}$ such that

[TABLE]

we have

[TABLE]

The condition (7h) is only required when the scatter density is a function of momentum transfer. For some modalities, scatter density varies with respect to scatter angle according to a known function (e.g. Klein–Nishina). In such cases the data can be corrected and momentum transfer removed from (7b).

For clarification, we will refer to the right-hand side of (7i) as the signed broken ray transform (SBRT) due to the sign change between CBTs. This is equivalent to the signed V-line transform [8]. Some authors have reserved their definition of the BRT for this later expression [2]. While either definition of the BRT assumes a linear combination of two CBTs sharing a common vertex, the distinction is important for inversion. Also, while positive images yield positive BRT data, SBRT data may be negative.

For tomographic imaging applications it is common to index the data according to the source and detector locations. In contrast our indexing is somewhat unconventional. In the context of the BRT, Katsevich and Krylov were the first to demonstrate the benefits of indexing the data by the scatter location [2]. Under this indexing schema, both the CBT and BRT are linear and shift-invariant (LSI). Linear systems analysis is therefore applicable to the CBT and BRT. Their relationship is depicted in Figure 1b. This is a central theme to our contribution as we are the first to consider the two-dimensional Fourier transform of the BRT. This perspective has benefits which we will demonstrate in subsequent sections.

Our focus is limited to 2D single-scatter imaging problems where scatter events are observed throughout the media of interest. This distinction is important because the terms broken ray transform and v-line transform have been used to describe a number of related problems. We distinguish BRT problems integrating over multiple reflections [17, 18] or integrating over multiple vertices [19, 20]. Some constrain the vertex locations along the perimeter of the measurement geometry [17, 18, 21]. This is generally motivated by the use of Compton cameras. In three dimensions this results in the cone transform [22], which we distinguish from the cone beam transform (4)[16] applicable to our measurement geometry.

The first analytic inversion formula for the BRT is due to Florescu et al. [4]. The global inversion formula requires only two scatter angles to recover the attenuation image in the presence of spatially varying scatter density. The inversion technique can be summarized as a three-step process. First, obtain the one-dimensional Fourier transform of the data. For the second step, each frequency is considered independently. Solve the resulting complex, one-dimensional, bounded differential equation. Third, obtain the inverse one-dimensional Fourier transform across the solutions. This yields an exact reconstruction of images with bounded support. The coordinate used to index the data in the original derivation were not linear-shift invariant. This was later derived using data indexed by the scatter location and generalized for higher dimensions [5].

The number of available scatter angles is a discriminating factor in selecting a BRT inversion strategy. A local inversion formula was discovered by Katsevich and Krylov requiring 3 unique scatter angles [2]. In contrast with prior results, [4], their reconstructions demonstrated significant reduction in artifacts. This was later generalized for additional scatter angles and source locations [6]. While the attenuation map can be recovered locally, the recovery of the scatter density image still requires global reconstruction of the attenuation image. This, and the requirement (7h) for coherent scatter imaging, motivates our interest in 2D BRT inversion techniques using only two scatter angles.

The initial results by Florescu et al. contained significant artifacts even for trivial phantoms [4]. These artifacts were broadly attributed to the nonlocal effects of integration. Artifacts in initial results exhibited striations at three distinct angles. Two of these angles are associated with the incident and scatter directions ($\theta_{i}$ and $\theta_{j}$). However, this did not directly address the third direction. This was later explored using micro-local analysis [7]. Sherson was the first to recognize $\theta_{i}+\theta_{j}$ as the direction of integration required for inversion.

Most recently a new inversion technique was developed by Ambartsoumian and Jebelli [8]. They used linear-shift invariant indexing of the data but did not employ the Fourier transform. They thoroughly and eloquently derive a new inversion technique by extending the Fundamental Theorem of Calculus to higher dimensions under a linear change of variables. They consider V-Line transformed (VLT) data defined by the linear combination CBTs along multiple directions $\{\theta_{i}\}$. Integrating VLT data along the direction $\sum_{i}\theta_{i}$ yields the integral of the image over the faceted cone defined by $\{\theta_{i}\}$ and the common vertex $x$. This unbounded volume can be reduced to a parallelepiped by linearly combining samples of this integral. Weighting the results, one obtains a reconstruction of the image averaged over this volume. This leads to a wonderfully concise inversion formula. In two dimensions they replace differentiation along the directions $\theta_{i}$ and $\theta_{j}$ with sample differences. This leaves only integration along the direction ${\theta_{i}+\theta_{j}}$. The consequence is potential blurring over the resulting parallelogram. This can be arbitrarily small for noise-free environments with high resolution data. For noisy data, the size of the parallelogram must be larger which effects blurring in the reconstruction. Additionally, artifacts appear along the direction of integration ($\theta_{i}+\theta_{j}$) [8].

In the following we take a fresh look at the BRT as a linear shift-invariant operator. A linear systems perspective provides new insights on the transforms and tools for contrasting prior inversion formulas. More specifically, we demonstrate images with bounded support do not guarantee data with bounded support. We are the first to consider the minimum data required for reconstruction and techniques for bounding support of the data. The two-dimensional Fourier transform of the BRT operator is especially useful for contrasting inversion techniques. It exhibits zeros along direction of integration ($\theta_{i}+\theta_{j}$) [8]. The ensuing ambiguity can be resolved analytically with boundary conditions on the reconstructed image (e.g. 0 after subtracting the background level). However, this does not address numerical sensitivity. The poles in the forward operator, along the directions $\theta_{i}$ and $\theta_{j}$, also present numerical challenges. We contrast recent work [8] against prior inversion techniques [4, 5, 7] as different strategies for addressing the poles in the forward operator. To mitigate numerical issues associated with the forward operator, we advocate Tikhonov regularization in obtaining reconstructed images. This can be implemented efficiently in the Fourier domain for arbitrary angles $\theta_{i}$, $\theta_{j}$ due to rotational invariance of the two-dimensional Fourier transform.

The remainder of this document is organized as follows. In Section 2 we provide new analysis of the BRT from a linear systems perspective. In Section 3 we present new numeric algorithms demonstrating some benefits of the preceding analysis. In Section 4 we present results of numerical simulations. Finally, some conclusions are given in Section 5.

2 New Analysis of the Broken Ray Transform

To exploit the benefits of linear system analysis, we first derive the two-dimensional Fourier transform of the BRT data. Since the BRT is LSI, we expect the result to have a specific form. We can decompose the transform of the data into the product of two terms: the transform of the image and the transform of the system function. From the transform of the system function, several insights are directly available.

We consider an absolutely integrable image with bounded support. We define a closed, bounded, convex set $C\subset\mathbb{R}^{2}$, which we use to window the image

[TABLE]

We use $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\mu}}}_{C}(w)$ to represent the two-dimensional Fourier transform of the image. Since the BRT is simply a linear combination of CBTs, we first define the two-dimensional Fourier transform of the CBT data

[TABLE]

The details of this derivation are in A. The two-dimensional Fourier transform of the BRT data is therefore

[TABLE]

Indeed, this result can be decomposed into the product of two terms. The bracketed term represents the two-dimensional Fourier transform of the BRT system function. For convenience, we will frequently reference a portion of this term

[TABLE]

We emphasize $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{h}}}_{i,j}(w)$ is not the transform of the BRT forward operator as the delta functions have been excluded.

The expression (7klc) highlights some challenges with BRT inversion. We observe singularities at $w\cdot\theta_{i}=0$ and $w\cdot\theta_{j}=0$. At these frequencies, finite $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\mu}}}_{C}(w)$ does not guarantee finite $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{g}}}_{i,j}(w)$. As a consequence, the data may have unbounded support.

Additionally (7klc) demonstrates zeros in the forward operator. We define the set $\Theta_{i,j}\subset\mathbb{R}^{2}$ as

[TABLE]

For all $w\in\Theta_{i,j}$, we have $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{g}}}_{i,j}(w)=0$ for all $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\mu}}}_{C}(w)$. In this way the BRT has a non-trivial nullspace. The zeros are limited to a line, and so the nullspace does not include images with bounded support. This does not preclude exact analytic reconstruction of images with bounded support. However, this is problematic for numeric reconstruction. We have arrived at these observations from a linear systems perspective. Similar observations were previously made applying microlocal analysis to the BRT [7].

A Fourier representation of the image is found multiplying both sides of (7klc) by the inverse of (7klm)

[TABLE]

Justification for removing the delta functions is given in B. Using (7klo) alone, we cannot recover $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\mu}}}_{C}(w)$ for $w\in\Theta_{i,j}$. According to (7klc), $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{g}}}_{i,j}(w)=0$, for all $w\in\Theta_{i,j}$, which leaves (7klo) indeterminate. This can be resolved imposing boundary conditions on $\mu_{C}(x)$ or, equivalently, continuity of $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{g}}}_{i,j}(w)$  (i.e. applying L’Hôpital’s rule).

In the parlance of linear systems analysis, this inversion formula comprises two lines of zeros and one line of pole. The zeros are associated with a directional derivatives, and the poles are associated with integration. Inverting this process leads to the reconstruction formulas

[TABLE]

A detailed derivation is given in B. This is a generalization of previous inversion formulas derived by other means [4, 5, 7, 8]. Sherson was the first to recognize the symmetry in this expression[7] which is useful for numeric reconstructions from noisy data. For finite data the integrals are applied over different lengths and different noise realizations. The two results can be combined to minimize variance in the reconstruction.

The delta functions in (7klc) present some obvious challenges. This implies images with bounded support do not guarantee data with bounded support. This begs the question: what extent of data is necessary for reconstruction? We address this in Section 2.1. Additionally, these present numerical challenges when computing the Fourier transform from sampled data. We will address this in Section 2.2. Combining these results we obtain Fourier-based inversion formulas in Section 2.3.

2.1 Complete Representation of Data With Infinite Support

It is helpful to distinguish segments of the boundary of $C$ with respect to the orthogonal basis $\theta$, $\theta^{\perp}$. For this we define the scalar values

[TABLE]

Additionally we define the auxiliary functions ${u^{-}_{\theta},u^{+}_{\theta}\mathrel{\mathop{\mathchar 58\relax}}\left[v^{-}_{\theta},v^{+}_{\theta}\right]\rightarrow\mathbb{R}}$

[TABLE]

We define non-overlapping line segments along the boundary of $C$ as functions, ${f^{-},f^{+}\mathrel{\mathop{\mathchar 58\relax}}\left[v^{-}_{\theta},v^{+}_{\theta}\right]\rightarrow C}$

[TABLE]

Using these functions, we define the mutually exclusive regions ${C^{-}_{\theta},C^{+}_{\theta},V^{-}_{\theta},V^{+}_{\theta}\subset\mathbb{R}^{2}}$:

[TABLE]

These definitions are illustrated in Figure 2a.

With definitions in place, we make some observations regarding the support of the CBT data. We state them as three theorems. First, we limit support of the data. The CBT data are zero for all $x$ outside the support of the image, $C$, and the shadow region $C^{-}_{\theta}$.

Theorem 1**.**

$(B\mu_{C})(x,\theta)=0$* for all $x\in\overline{C\cup C^{-}_{\theta}}$*

Proof.

The set $\overline{C\cup C^{-}_{\theta}}$ can be partitioned into three regions $C^{+}_{\theta}$, $V^{+}_{\theta}$, and $V^{-}_{\theta}$, where

[TABLE]

Since $C^{+}_{\theta}\cap C=\varnothing$, we have $\mu_{C}(x)=0$ for all $x\in C^{+}_{\theta}$. Additionally, if $x\in C^{+}_{\theta}$, then according to the definition (7klpy), ${x+s\theta\in C^{+}_{\theta}}$ for all $s\geq 0$. Therefore ${(B\mu_{C})(x,\theta)=0}$ for all $x\in C^{+}_{\theta}$. The same is true for $V^{-}_{\theta}$ and $V^{+}_{\theta}$. ∎

Next we observe the CBT, over the region $C_{\theta}^{-}$, is constant along the direction $\theta$. The values are determined by the Radon transform at $\theta$.

Theorem 2**.**

$(B\mu_{C})(x,\theta)=(R\mu_{C})(x\cdot\theta^{\perp},\theta)$* for all $x\in C^{-}_{\theta}$*

Proof.

For all $x\in C^{-}_{\theta}$, we can extend the integral of the CBT

[TABLE]

where the first term is 0 due to the bounded support of $\mu_{C}$. Combining these integrals and expanding $x$ along the orthogonal basis vectors $\theta$ and $\theta^{\perp}$, we have

[TABLE]

∎

Finally, the Radon transform, for fixed direction $\theta$, is given by the CBT along the boundary of $C$.

Theorem 3**.**

${(R\mu_{C})(v,\theta)=(B\mu_{C})(f^{-}_{\theta}(v),\theta)}$* for all ${v\in\left[v^{-}_{\theta},v^{+}_{\theta}\right]}$ and ${(R\mu_{C})(v,\theta)=0}$ for all ${v\notin\left[v^{-}_{\theta},v^{+}_{\theta}\right]}$.*

Proof.

For ${v\in\left[v^{-}_{\theta},v^{+}_{\theta}\right]}$, we can expand

[TABLE]

For $v<v^{-}$, we have ${\mu_{C}(v\theta^{\perp}+t\theta)=0}$, since ${v\theta^{\perp}+t\theta\in V^{-}_{\theta}}$ for all $t$. Therefore $(R\mu_{C})(v,\theta)=0$ for $v<v^{-}_{\theta}$. The same can be shown for $v>v^{+}_{\theta}$. ∎

Corollary 3.1**.**

$(B\mu_{C})(x,\theta)=(R\mu_{C})(x\cdot\theta^{\perp},\theta)$* for all ${x\in V^{-}_{\theta}\cup C^{-}_{\theta}\cup V^{+}_{\theta}}$*

Proof.

Theorem 2 demonstrates equality for $x\in C^{-}_{\theta}$. For ${x\in V^{-}_{\theta}\cup V^{+}_{\theta}}$, we have $(B\mu_{C})(x,\theta)=0$ according to Theorem 1. For ${x\in V^{-}_{\theta}}$, we have $x\cdot\theta^{\perp}<v^{-}_{\theta}$. Therefore, ${(R\mu_{C})(x\cdot\theta,\theta)=0}$ according to Theorem 3. Similarly, ${(R\mu_{C})(x\cdot\theta,\theta)=0}$ for ${x\in V^{+}_{\theta}}$. ∎

Theorem 2 demonstrates images with bounded support do not guarantee CBT data with bounded support since $C^{-}_{\theta}$ is unbounded. This is problematic for discrete Fourier analysis. However, data outside the support of the original image is redundant. If $(B\mu_{C})(x,\theta)$ is known for all $x\in C$ including its boundary, $(R\mu_{C})(v,\theta)$ is available. Combining Theorem 1 and Theorem 2, the CBT is then known for all $x\in\mathbb{R}^{2}$. This is significant as there may be problems for which data are not available outside the support of the original image. This demonstrates samples along the boundary, or alternatively direct-path (ballistic) measurements, are sufficient. Once this minimum extent of data are available, CBT data can be extended arbitrarily.

For our problems of interest CBT data are not available directly. The BRT is a linear combination of two CBTs sharing a common vertex. Similar to CBT data, bounded support of the image does not guarantee bounded support of the BRT data. The previous analysis of CBT data informs the sampling requirements on BRT data. Using the definition of the BRT (5), we distinguish two directions $\theta_{i}\neq\theta_{j}$. In addition to knowing $(G\mu_{C})(x,\theta_{i},\theta_{j})$ for all $x\in C$, we additionally require the Radon transform in two direction: $(R\mu_{C})(v,\theta_{i})$, and $(R\mu_{C})(v,\theta_{j})$. The complication lies in the partitions of the BRT data. BRT data requires additional partitions which may overlap. Resolving the Radon transform, with respect to two directions, from the BRT data is more challenging.

Following the previous work, our definitions for $C$ and $\mu_{C}$ need not change. However, we use $i$ and $j$ to distinguish the directions the subscripts of the definitions (7klpq)-(7klpaa). These indices are used only in subscripts to avoid confusion with the imaginary unit $i\mathrel{\mathop{\mathchar 58\relax}}=\sqrt{-1}$. Depending on $C$, $\theta_{i}$, and $\theta_{j}$, the set $C_{i}^{-}\cap C_{j}^{-}$ may be nontrivial. The BRT data can be partitioned

[TABLE]

These regions are depicted in Figure 2b.

In contrast to the CBT, we must distinguish $(R\mu_{C})(v,\theta_{i})$ from $(R\mu_{C})(v,\theta_{j})$. Over $C$ alone, they may not be directly available. We consider two scenarios. First, for some regions $C$ and scatter angles $\theta_{i},\theta_{j}$, the set $C_{i}^{-}\cap C_{j}^{-}$ is empty. For example, this is true for rectangular $C$, when $\theta_{i}$ is parallel to a boundary of $C$, and $\theta_{i}\cdot\theta_{j}\leq 0$. In such cases, $(R\mu_{C})(v,\theta_{i})$ from $(R\mu_{C})(v,\theta_{j})$ can be distinguished along the boundary of $C$. As a second scenario, forward scatter (ballistic) measurements at the two angles can be used to measure the Radon transforms directly. This would require new measurements. However, this may be useful for some modalities if measurements over the boundary of $C$ are not available.

The notation introduced in this section is also useful for simplifying the assumed support of the image. Due to shift-invariance of the BRT, we can assume the image is centered about the origin without loss of generality.

Definition 2.1**.**

Let $C$ represent a closed and bounded region in $\mathbb{R}^{2}$, and let $\theta_{i}$ and $\theta_{j}$ represent unique directions such that $|\theta_{i}\cdot\theta_{j}|<1$. We define $v_{i}^{-}$, $v_{i}^{+}$, $v_{j}^{-}$, $v_{j}^{+}$ using (7klpq) and (7klpr). Then, $C$ is centered with respect to $\theta_{i}$ and $\theta_{j}$ when both ${v_{i}^{+}=-v_{i}^{-}}$ and ${v_{j}^{+}=-v_{j}^{-}}$.

Parallelograms are an important geometric shape in the context of the BRT. This was first recognized by Ambartsoumian and Jebelli [8]. It is often convenient to extend $C$ to the circumscribed parallelogram.

Definition 2.2**.**

Let $C$ represent a closed and bounded region centered with respect to $\theta_{i}$ and $\theta_{j}$. The circumscribed parallelogram, with edges parallel to $\theta_{i}$ and $\theta_{j}$, is given by

[TABLE]

where $v_{i}$, $v_{j}$ are defined according to (7klpq)-(7klps). Clearly, $C\subseteq P\subset\mathbb{R}^{2}.$.

In (7klpag), $P$ is expressed in terms of the orthogonal distance between parallel sides. Alternatively, we obtain the edge lengths $\alpha_{i}$ and $\alpha_{j}$ for the edges parallel to $\theta_{i}$ and $\theta_{j}$, respectively

[TABLE]

We can equivalently express $P$ in terms of the edge lengths

[TABLE]

Related to $P$, we define a parallelogram indicator function in C and derive its two-dimensional Fourier transform. The results will be referenced frequently in subsequent sections.

2.2 Filtering Unbounded Support of the Data

When the Fourier transform must be determined numerically, unbounded support of the BRT is problematic. Simply truncating BRT data effects blurring in the frequency domain. This corrupts the spectral representation and invalidates the previous Fourier reconstruction methods. Alternatively, we consider convolving data in the spatial domain. We define a generalized point spread function (PSF) such that the shifted copies of the data combine destructively outside a bounded region of support.

We consider the PSF

[TABLE]

Here $a_{i},a_{j}>0$, determine the shift lengths. The expression (7klpak) has the Fourier transform

[TABLE]

To reduce the number of variables defined we introduce new notation to distinguish signals, which support expansion using the PSF function (7klpak). We define

[TABLE]

The same superscript $m$ will be subsequently applied to continuous signals in the spatial domain, and sampled signals. Plugging (7klc) and (7klpak) into (7klpam) we have

[TABLE]

The inverse two-dimensional Fourier transform of this expression involves integration over $w$. Due to the sampling property of the delta function, and since $\sin(0)=0$, the final two bracketed terms vanish under integration. By the uniqueness of the Fourier transform, we have

[TABLE]

In (7klpapb) we make use of both (7klpan), and (7klm). Since the BRT is LSI, this result is expected. Filtering the input to an LSI system is equivalent to filtering the output. The significance is that the delta functions vanish when we filter the data using (7klpak).

We obtain another useful form by expanding (7klpapa) using (7klpal). We reappropriate the denominator of $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{h}}}_{i,j}(w)$ to find

[TABLE]

The product of $\operatorname{sinc}$ functions in (7klpapaq) is associated with a parallelogram window function as demonstrated in C. This motivates the definition

[TABLE]

where $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{p}}}_{i,j}\left(w;a_{i},a_{j}\right)$ is defined according to (7klpapcy). The scaling is motivated by (7klpapcp). Using (7klpapar) in (7klpapaq), we have

[TABLE]

Taking the inverse two-dimensional Fourier transform of (7klpapas) we find

[TABLE]

Here the first term represents the directional derivative in the direction $\theta_{i}+\theta_{j}$. This is clearly not a unit vector. In this form we observe $g^{m}_{i,j}(x;a_{i},a_{j})$ has bounded support.

Theorem 4**.**

For an absolutely integrable image with bounded support, filtering the BRT data with the PSF (7klpak) bounds support of the data for all $a_{i},\,a_{j}\in\left(0,\infty\right)$. Additionally, the data are finite everywhere.

Proof.

Without loss of generality, we assume the support of the image $\mu_{C}(x)$ is bounded by the circumscribed parallelogram, $P$, according to Definition 2.2. We first observe $\mu_{i,j}^{p}\left(x;a_{i},a_{j}\right)$ has bounded support. We define

[TABLE]

The indicator function $p_{i,j}\left(x;a_{i},a_{j}\right)$, defined by (7klpapco), has bounded support over a parallelogram similar to $P$ in (7klpag). Taking the convolution of two functions defined over similar parallelograms, the support of the result is also bounded by a similar parallelogram. This limits support of $f(x)$ to a parallelogram with sides parallel to $\theta_{i}$ and $\theta_{j}$ with perpendicular distances $v_{i}+b_{i}$ and $v_{j}+b_{j}$, respectively. The variables $b_{i}$ and $b_{j}$ are related to $a_{j}$ and $a_{i}$ according to (7klpapcl) and (7klpapcm), respectively.

Using (7klpapaub) in (7klpapat), we have

[TABLE]

Outside the the region of support of $f(x)$, its directional derivative is also zero. Therefore $g^{m}_{i,j}(x;a_{i},a_{j})$ has bounded support.

To show $g^{m}_{i,j}(x;a_{i},a_{j})$ is finite everywhere, we consider

[TABLE]

This is finite due to the assumption $\mu_{C}(x)$ is integrable. ∎

2.3 Image Reconstruction from BRT Data with Bounded Support

Bounded BRT data facilitates numeric inversion in the frequency domain. We consider two inversion strategies. The size of the available spreading parameters $a_{i}$ and $a_{j}$ in (7klpak) plays an important role in selecting an inversion strategy. In both cases we reconstruct a version of the desired image subject to convolution. However, the PSFs associated with the reconstructed images are different.

Multiplying both sides of (7klpapb) by the inverse of (7klm) we have the relationship

[TABLE]

This is similar to (7klo). However, the reconstruction is subject to multiplication with the PSF (7klpal). Analytically, we can recover $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\mu}}}^{m}_{i,j}\left(w;a_{i},a_{j}\right)$ from $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{g}}}^{m}_{i,j}(w;a_{i},a_{j})$ using (7klpapax) and continuity assumptions or, equivalently, boundary conditions on $\mu^{m}_{C}(x)$.

We find a representation of the left hand side of (7klpapax) in the spatial domain by taking the inverse two-dimensional Fourier transform of (7klpan)

[TABLE]

For small $a_{i}$ and $a_{j}$, the image copies will overlap. As $a_{i}$ and $a_{j}$ increase we can reconstruct $\mu_{C}\left(x\right)$ from segments without overlap.

Theorem 5**.**

An image with with bounded support, $\mu_{C}(x)$, can be recovered from filtered BRT data $g^{m}_{i,j}(x;a_{i},a_{j})$ when ${a_{i}>v_{j}/\left|2\det\left(\theta_{i},\theta_{j}\right)\right|}$ and ${a_{j}>v_{i}/\left|2\det\left(\theta_{i},\theta_{j}\right)\right|}$ for $v_{i}$, $v_{j}$ defined according to (7klps).

Proof.

A portion of the image $\mu_{C}(x)$, without overlap, is associated with each shifted copy in (7klpapay). When the shifts are sufficiently large, the partial images can be combined to reconstruct the original image. To demonstrate this it is useful to first extend $C$ to the circumscribed parallelogram $P$ in (7klpaj). The edge lengths of this geometric region bound the minimum shift lenghts for image recovery.

To emphasize $P$ as the assumed region of support, we use $\mu_{P}(x)$. Since $C\subseteq P\subset\mathbb{R}^{2}$, we have $\mu_{P}(x)=\mu_{C}(x)$, for all $x\in\mathbb{R}^{2}$. For $a_{i}>\alpha_{i}/2$ and $a_{j}>\alpha_{j}/2$ we can recover $\mu_{P}(x)$ from $\mu^{m}_{i,j}(x)$ using

[TABLE]

Each case can be expanded as a series of four terms using (7klpapay). However, three of these terms are zero due to the support of $P$ in (7klpaj). Combining the four cases we recover $\mu_{P}(x)$, and therefore $\mu_{C}(x)$, for all $x\in\mathbb{R}^{2}$. Expanding $a_{i}>\alpha_{i}/2$ and $a_{j}>\alpha_{j}/2$ using (7klpah), (7klpai) we obtain the boundary in the stated form. ∎

In general this approach requires large data sets to obtain $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{g}}}^{m}_{i,j}(w;a_{i},a_{j})$. For many cases we can extend BRT data using the techniques in Section 2.1. However, when we are limited to small $a_{i}$ and $a_{j}$, another approach is necessary.

Alternatively, we can simply recover $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\mu}}}_{i,j}^{p}\left(w;a_{i},a_{j}\right)$. From (7klpapas), we have

[TABLE]

Taking the inverse two-dimensional Fourier transform of (7klpapba) we have

[TABLE]

This is equivalent to the inversion formula of Ambartsoumian and Jebelli [8]. For sampled data this formula can be implemented easily whenever the direction of integration is aligned with a sampling axis. For other cases, the frequency domain representation (7klpapba) is useful.

We emphasize $\mu_{i,j}^{p}\left(x;a_{i},a_{j}\right)\neq\mu_{C}\left(x\right)$. Taking the inverse two-dimensional Fourier transform of (7klpapar), we have

[TABLE]

This demonstrates the recovered image as a blurring of the original image with a parallelogram window function. For high resolution, noise-free, data the size of this window can be made arbitrarily small. The recovery (7klpapbb) is only approaches $\mu_{C}(x)$ in a limiting sense [8].

3 Numeric Algorithms

In application we must reconstruct images from sampled data. Our analysis of the BRT from a linear systems perspective extends easily to sampled data. We demonstrate this with two new algorithms. First we provide an algorithm for extending CBT data. This is motivated by the work in Section 2.1. For a broad class of problems this can be applied to BRT data and therefore facilitates use of the filtering methods of Section 2.2. For brevity, we do not address numeric implementation of the filtering methods. However, D is useful in this context. Second, we present a numeric inversion algorithm for bounded BRT data. Leveraging the rotational invariance of the two-dimensional Fourier transform, the directions $\theta_{i},\theta_{j}$ are unconstrained in our algorithm. We also include regularization to address poor conditioning of the forward operator.

3.1 Extending Cone Beam Transform Data

We consider CBT data sampled uniformly over a rectangular region. For consistency with previous definitions, we expand $x$ along two scalar axes $x=(t,y)$. For the two axes we use subscripts to distinguish the number of samples $N_{t}$, $N_{y}$ and the sample spacing $\Delta_{t}$, $\Delta_{y}$. We collect the available data in the $N_{y}\times N_{x}$ matrix $B$. The elements represent samples of the CBT data

[TABLE]

for $n\in\{1,\ldots,N_{y}\}$, and $m\in\{1,\ldots,N_{t}\}$. We expand the direction $\theta=(\cos\xi,\sin\xi)$. The spatial location associated with sample $B_{1,1}$ is $x_{B}=(t_{B},y_{B})$. In this configuration the $y$ coordinate increases with the row index $n$, and the $t$ coordinate increases with the column index $m$. It is not necessary to distinguish the terms $\Delta_{t}$, $\Delta_{y}$, and $\xi$ for most of the computations related to sampled CBT data. For convenience we define

[TABLE]

which is a sufficient input for algorithms on uniformly sampled data.

Extending CBT data in the direction $\theta$ is trivial. For this we need only consider ${x\in V^{-}_{\theta}\cup C^{+}_{\theta}\cup V^{+}_{\theta}}$, where $(B\mu_{C})(x,\theta)=0$ according to Theorem 1. Zero padding is sufficient.

Extending the data in the direction $-\theta$ is nontrivial. For simplicity we first consider only $\xi\in(0,\pi/2)$. Figure 3 illustrates the problem of extending the data, $B$, into the quadrants $Q2$, $Q3$, and $Q4$. The Radon transform serves as a proxy for extending the data according to Corollary 3.1. We assume the first row and column comprise no samples interior to $C$ such that these data are samples of ${(R\mu_{C})(x\cdot\theta^{\perp},\theta)}$. We can then extend the data using $(B\mu_{C})(x,\theta)=(R\mu_{C})(x\cdot\theta^{\perp},\theta)$. A brute-force approach would be to resample the Radon transform for each new data point. A computationally efficient approach is to extend the CBT data by shifting samples along the boundaries. This process is detailed in Algorithm 1. For $\xi\notin(0,\pi/2)$, we can still use Algorithm 1 by suitably flipping the inputs and outputs.

We can adapt this process to an important class of BRT problems. We consider the incident direction aligned with the $t$-axis and $\theta_{i}\cdot\theta_{j}<0$. Specifically, we use $\theta_{i}=(-1,0)$ and expand $\theta_{j}=(\cos\xi,\sin\xi)$. We construct the BRT data matrix $G\in\mathbb{R}^{N_{y}\times N_{t}}$ with elements

[TABLE]

The previous definitions for $\Delta_{t}$, $\Delta_{y}$, and $\lambda$ remain applicable.

Extending BRT data requires knowledge of both $(R\mu_{C})(v,\theta_{i})$ and $(R\mu_{C})(v,\theta_{j})$. We assume BRT data are sampled beyond the support of the image, such that no boundary samples of $G$ correspond to points within $C$. For $\theta_{i}=(-1,0)$, and ${|\xi|<\pi/2}$ this implies ${[G]_{1,N_{t}}=[G]_{N_{y},N_{t}}=0}$. In this case $(R\mu_{C})(v,\theta_{i})$ can be recovered from the last column of $G$. We can extend BRT data in the direction $-\theta_{i}$ simply by repeating the last column. For $\xi>0$, the last row (maximum $y$) of $G$ is then zero. The function $(R\mu_{C})(v,\theta_{j})$ can be recovered from the first column of $G$ and the first row. The BRT data can be extended in the direction $-\theta_{j}$ using Algorithm 1. Alternatively, for $\xi<0$, $(R\mu_{C})(v,\theta_{j})$ can be recovered from the first column of $G$ and the last row. The BRT data can still be extended in the direction $-\theta_{j}$ using Algorithm 1. However, the inputs and outputs must be flipped accordingly.

3.2 Inversion of BRT Data with Regularization

Filtering ensures bounded support of $g^{m}_{i,j}(x;a_{i},a_{j})$. However, recovery of $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\mu}}}^{m}_{i,j}(w;a_{i},a_{j})$ is still ill-posed due to conditioning of $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{h}}}_{i,j}(w)$. For this we use Tikhonov regularization which can be applied sample-wise in the frequency domain.

We restate (7klm) as an expression of scalar values by expanding $w=(w_{t},w_{y})$, $\theta_{i}=(\cos\xi_{i},\sin\xi_{i})$, and $\theta_{j}=(\cos\xi_{j},\sin\xi_{j})$

[TABLE]

Notice this expression is commutative with respect to $\xi_{i}$ and $\xi_{j}$. We define the system matrix $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{H}}}$ by sampling (7klpapbg) uniformly

[TABLE]

The discrete analog of (7klpapb) is then

[TABLE]

Here we have used $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{G}}}^{m}$ to represent the two-dimensional discrete Fourier transform of $G^{m}$, the filtered analog of (7klpapbf). We use $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\Psi}}}^{m}$ to represent samples of $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\mu}}}^{m}_{i,j}(w)$. The symbol $\odot$ represents element-wise multiplication.

Zeros in the denominator of (7klpapbg) are problematic for numeric analysis. We define the auxiliary function

[TABLE]

However, filtering ensures $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{G}}}^{m}$ and $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\Psi}}}^{m}$ are also zero when ${d(w_{t},w_{y})=0}$. Zeros in the numerator of (7klpapbg) also affect conditioning of the problem. Tikhonov regularization provides a generic mitigation strategy. Putting this together, we approximate the element-wise inverse of $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{H}}}$

[TABLE]

where ∗ indicates complex conjugation, and $\epsilon$ is the smoothing parameter. This yields the estimate

[TABLE]

Applying the 2D inverse discrete Fourier transform to the result we obtain a reconstruction of the filtered attenuation image. This process is described in Algorithm 2. The smoothing parameter $\epsilon$, in (7klpapbk), can be adjusted for measurement noise and numerical errors.

Tikhonov regularization is generic and does not impose boundary conditions. For arbitrary angles $\xi_{i}$, and $\xi_{j}$, few samples of $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\Psi}}}^{m}$ lie in the nullspace of the forward operator and it is sufficient to zero the results at these samples. Otherwise, it may be necessary to impose boundary conditions. For example, to ensure $\Psi^{m}$ is zero along the $t=0$ and $y=0$ boundary, all columns and all rows of $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\Psi}}}^{m}$ must sum to 0.

4 Numerical Simulations

We provide results of numerical simulations to demonstrate the utility of this analysis. We use the modified Shepp-Logan phantom [23, 24] in most of our simulations as depicted in Figure 4. This phantom is reasonably challenging and the BRT data can be determined analytically. For Figure 4 we sample the image and data space uniformly in $y$ and $t$. For $y$ we use $N_{y}=600$ sampling over $[-1,1]$. For $t$ we use $N_{t}=400$ sampling over $[-0.75,0.75]$. This effects different sampling rates in $t$ and $y$. Limiting the extent of available BRT data in this way truncates the data both in $y$ and $t$ as shown in Figure 4b and Figure 4c

We first demonstrate filtering bounds support of the data. Results are shown for both the BRT and SBRT in Figure 5. In this case the filtered image and filtered data were all obtained analytically and then sampled.

Filtering can also be applied to sampled BRT data directly. For sampled data this effects small errors which we quantify against the reference data of Figure 5. Results are shown in Figure 6. Artifacts are observed at scatter points for which resulting rays are tangent to large transitions in the image. This is a consequence of sampling. For both BRT and SBRT filtering the peak absolute error is less than 5% the peak image value.

Further analysis of $h(w)$ provides insights on BRT inversion. We can express (7klm) in polar coordinates with the change of variables

[TABLE]

such that

[TABLE]

We make a few observations. First, $\rho$ in the denominator of (7klpapbn) implies the BRT attenuates high frequency content. Reconstruction will be sensitive to noise at high frequencies. Second, there are singularities at ${\phi=\xi_{i}\pm\pi/2}$ and ${\phi=\xi_{j}\pm\pi/2}$. Filtering ensures the image and data are zero at these frequencies. Finally, (7klpapbn) is zero at ${\phi=\frac{1}{2}\left(\xi_{i}+\xi_{j}\right)\pm\pi/2}$. These zeros do not appear in the CBT, but arise in the combination of two CBTs.

The matrix $K$ plays a critical role in BRT reconstruction (7klpapbl). This incorporates changes to $h(w)$ due to $\xi_{j}$, and the regularization term $\epsilon$. Changes to $|K|$ with respect to these terms is shown in Figure 7. Here we fix ${\xi_{i}=\pi}$ without loss of generality. The lines indicating strong attenuation are due to singularities of (7klpapbn) at ${\phi=\pm\pi/2}$, and ${\phi=\xi_{j}\pm\pi/2}$. Zeros in (7klpapbn) effect large amplitudes in $|K|$ along ${\phi=\xi_{j}/2}$. However, this amplitude is curtailed through regularization as $\epsilon$ increases. Without regularization, we would expect reconstruction artifacts along this spectral line.

The original global BRT inversion formula is due to Florescu et al. [4]. We will refer to this as the FMS formula. Specifically contrasting with our algorithm, we analyze the same square phantom in Figure 8. The original work assumed data available over an infinite strip with no additional insights on limiting the data. The data of Figure 8b violates this assumption. Directly applying the FMS formula to this data yields poor results as shown in Figure 8c. However, we can simulate additional data using Algorithm 1. Applying the FMS formula to the extended BRT data yields results consistent with those previously published [4]. In this way, Algorithm 1 can be used as a preprocessing step to reduce the extent of sampling required for reconstruction using the FMS formula. The direction of the artifacts is explained by the nullspace of the forward operator (7klm). Striations are observed in the direction ${\xi/2+\pi/2}$. Regularization further improves reconstruction as demonstrated in Figure 8f.

Inversion results for the Shepp-Logan phantom on noisy data are shown in Figure 9. The BRT data were obtained analytically, sampled, and corrupted with additive Gaussian noise. For small $\epsilon$, we see artifacts where the direction $\xi$ is tangent to high frequency edges of the image. This is a consequence of sampling errors and extending the BRT data. Additionally, edges perpendicular to the direction $\xi/2$ are not well resolved. This effects blurring along the direction $\xi/2$. Increasing $\epsilon$ increases the angular extent of blurring. The effect is reduced as $\xi$ increases.

5 Discussion

In conclusion we have demonstrated a new inversion algorithm for the BRT with improved performance with one scatter angle. Casting the BRT as a linear combination of CBTs provides insight on the minimum extent of sampling required for reconstruction and techniques to bound support of the data. Indexing the data by the scatter location, the BRT is an LSI operator. Analyzing the BRT as a linear operator, in the Fourier domain, yields a concise representation of its nullspace and highlights numerical sensitivities. This motivates the use of regularization which improves reconstruction from sampled data. Improving inversion of the BRT using one scatter angle supports extension to coherent scatter x-ray imaging problem which is angularly selective due to the momentum transfer.

• •

J. A. O’Sullivan was supported by NIH R01 CA 212638.

• •

This paper has been improved substantially by the reviewers’ comments.

Appendix A Derivation of the Fourier Transform of the CBT

We define the Fourier transform as a function of frequency to avoid scaling the inverse. For a one-dimensional function we define the one-dimensional Fourier transform and its inverse

[TABLE]

For two-dimensional functions we define the two-dimensional Fourier transform and its inverse

[TABLE]

In this form, we have

[TABLE]

where $\delta(x)$ and $u(x)$ represent the Dirac delta function and the unit step function, respectively.

For CBT data associated with a fixed direction, $\theta$, we define the two-dimensional Fourier transform

[TABLE]

In (7klpapbv) we changed the order of integration and substituted ${y=x+t\theta}$. In (7klpapbw) we substituted ${\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{\mu}}}_{C}(w)=\mathcal{F}^{2}\{\mu_{C}(x)\}}$. Finally, in (7klpapbx) we made use of (7klpapbs).

Appendix B BRT Inversion by Fourier Analysis

To invert the BRT, we start by multiplying both sides of (7klc) with the reciprocal of (7klm). Rearranging terms, we have

[TABLE]

The first two terms on the right hand side vanish under integration. We note the inverse two-dimensional Fourier transform

[TABLE]

Incorporating multiplicative functions does not change this result as long as they are finite for all $w\cdot\theta=0$. For the first term in (7klpapby), we expand ${w=s\theta_{i}+t\theta_{i}^{\perp}}$ as an orthonormal basis and set $s=0$. This leads to

[TABLE]

which is finite for all $t$ by our assumptions on $\mu_{C}(x)$. Applying a similar process for the second term, we find

[TABLE]

To address the third term in (7klpapby) we make use of the derivative and integral properties of the Fourier transform. Now we expand $x=s\theta+t\theta^{\perp}$ and consider the directional derivative

[TABLE]

Applying this to the inverse two-dimensional Fourier transform, we find

[TABLE]

We previously derived the integration property of the two-dimensional Fourier transform in the context of the CBT. From (7klpapbx), we have

[TABLE]

When $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{f}}}(w\cdot\theta)=0$ for all $w\in\mathbb{R}^{2}$, substituting $-\theta$ in (7klpapcg) we also find

[TABLE]

From (7klc), $\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{g}}}_{i,j}(0)$ is not guaranteed to be finite, much less zero. However, $\left(w\cdot\theta_{i}\right)\left(w\cdot\theta_{j}\right)\hat{\vphantom{\rule{1.0pt}{6.36115pt}}\smash{\hat{g}}}_{i,j}(w)=0$ for all $w\cdot\left(\theta_{i}+\theta_{j}\right)=0$. Putting this all together, we have equivalent expressions

[TABLE]

We emphasize equality only holds for images, $\mu_{C}(x)$, with bounded support. The assumption is necessary due to the nullspace of the forward operator.

Appendix C Two-Dimensional Fourier Transform of a Parallelogram

Parallelograms are often expressed in terms of the edge directions and edge lengths. We consider the directions $\theta_{i}$, $\theta_{j}$ and associated edge lengths $a_{i}$, $a_{j}$, respectively. The total area of the resulting parallelogram is ${a_{i}a_{j}|\det\left(\theta_{i},\theta_{j}\right)|}$. As an alternative to edge lengths, we also consider the orthogonal distance between parallel sides. We define $b_{i}$ and $b_{j}$ as the extent (height) of the parallelogram in the orthogonal directions $\theta_{i}^{\perp}$ and $\theta_{j}^{\perp}$, respectively. These distances are related to the edge lengths through the change of variables

[TABLE]

Additionally, we define the one-dimensional rectangular function

[TABLE]

We define the two-dimensional parallelogram indicator function, centered at $x=0$,

[TABLE]

Here $b_{i}$, $b_{j}$ are determined by $a_{j}$, $a_{i}$ using (7klpapcl) and (7klpapcm), respectively. The area of this function is equivalently

[TABLE]

To determine the two-dimensional Fourier transform of (7klpapco), we exploit the convolution property of the Fourier transform. We transform the two rectangular functions separately, then convolve the results in the frequency domain. The one-dimensional Fourier transform of (7klpapcn)

[TABLE]

Extending this to two dimension, we have the relation

[TABLE]

We derive the two-dimensional Fourier transform of (7klpapco)

[TABLE]

For (7klpapcv), we expand the integration variable in (7klpapcu) using the orthonormal basis, ${y=s\theta_{i}+t\theta_{i}^{\perp}}$, and integrate over $s$. Changing the variable of integration again effects a change in scaling in (7klpapcw). The second $\operatorname{sinc}$ function of (7klpapcw) comprises expansion of $\theta_{i}$ using the orthonormal basis $\theta_{j}$, $\theta_{j}^{\perp}$. Restoring $\theta_{i}$, we obtain (7klpapcx). Restoring $a_{j}$, $a_{i}$ using (7klpapcl), (7klpapcm), we obtain (7klpapcy).

Appendix D Non-Integer Shifts of Sampled Signals

Non-integer shifts of sampled signals requires interpolation. Fast implementations of the discrete Fourier transform (DFT) can be leveraged to perform this task quickly. For continuous signal $x(t)$, and uniform sample spacing $\Delta$, we define the $N$-length sampled signal

[TABLE]

The Fourier coefficients are given using the DFT

[TABLE]

for $m\in\{0,\ldots,N-1\}$. The corresponding inverse DFT is

[TABLE]

Using the generative property of the Fourier series representation[25], we can approximate

[TABLE]

Since $s$ is represented in samples, equation (7klpapde) is independent of sampling rate. This is particularly efficient when multiple shifted copies of the same signal are required. In such cases $y[m]$ need only be computed once. Additional savings are realized computing the $\operatorname{DFT}^{-1}$ in (7klpapde) for all signals at once. This process is described in Algorithm 3. Here we have included additional inputs indicating zero padding, $p$, and fill samples ${\bf f}$ to reduce aliasing.