Keywords: V-line Transforms, Tensor Tomography, Inversion Algorithms, Numerical Implementation
Mathematics subject classification 2020: 44A12, 44A30, 44A60, 47G10, 65R10, 65R32

1 Introduction↩︎

Optical tomography has emerged as an important and rapidly developing imaging modality in modern imaging science. It uses transmitted and scattered light to obtain information about an object’s internal structure and physical properties. In most imaging applications, the objective is to reconstruct spatially varying absorption and scattering coefficients from measurements collected at the object’s boundary. Such reconstruction problems naturally lead to inverse problems involving integral geometric transforms. Under suitable physical assumptions, it is commonly assumed that most photons undergo only a single scattering event while propagating through the object (see [3]–[5]). This single-scattering behavior naturally leads to the terminology broken-ray or V-line transform.

Over the past several years, significant attention has been devoted to imaging models based on generalized Radon transforms, particularly the V-line transform. These transforms generalize the classical Radon transform by integrating functions or tensor fields over broken-ray trajectories composed of two straight segments that intersect at a common vertex. Depending on the underlying physical model, the unknown quantity may be represented by a scalar function, vector field, symmetric \(2\)-tensor field, or, more generally, a higher-order symmetric tensor field. Consequently, the inversion and analysis of these V-line transforms are fundamental for the corresponding image reconstruction problems. V-line transforms naturally arise in several imaging modalities, including single-scattering optical tomography [3], [4], single-scattering X-ray tomography [6], [7], fluorescence imaging [8], Compton camera imaging [9], [10], Compton-scattering emission tomography with collimated detectors [11], [12], and gamma-ray transmission/emission imaging [13].

Due to their broad range of applications, V-line transforms have been extensively investigated from both theoretical and computational points of view; see, for example, [14]–[26] and the references therein. Previous studies on V-line transforms have primarily focused on two distinct geometric configurations. The first class consists of V-lines whose vertices lie on or outside the boundary of the imaging domain; see, for instance, [27], [28]. The second class consists of V-line trajectories whose vertices are located within the support of the unknown function. Such transforms naturally arise in single-scattering tomography and Compton camera imaging. The present work concerns this second class of V-line transforms, and we briefly discuss relevant results in the next paragraph.

In [29], [30], the authors investigated V-line transforms for vector fields and symmetric \(2\)-tensor fields in a setting where the vertex moves inside the disk while the V-line directions remain fixed. They derived exact kernel characterizations for these transforms and showed that vector fields and symmetric \(2\)-tensor fields can be uniquely reconstructed from appropriate combinations of V-line transform measurements and their corresponding moments. The numerical implementation, verification, and validation of these reconstruction algorithms were later developed in [31], [32].

A different geometry of V-lines whose vertex located inside the disk, described in Section 2, was introduced by [33]. In that work, the author derived theoretical and numerical reconstruction formulas for the V-line Radon transform of scalar functions under the support restriction that lies inside a disk. Subsequently, in [34], a complete reconstruction procedure was established based on a series-type inversion formula for the same transform. The numerical verification of these inversion methods was further studied in [35]. For a detailed account of the mathematical foundations, historical development, and applications of the generalized Radon transforms, we refer the reader to the recent monograph [36] by Gaik Ambartsoumian. The works [1], [2] extended the reconstruction results of [33], [34] from scalar functions to vector fields and higher-order symmetric tensor fields respectively. In these studies, inversion formulas were derived for the recovery of vector fields and symmetric tensor fields from appropriate collections of V-line transform data.

Motivated by all these developments, the present article makes several contributions. First, we establish an explicit characterization of the kernels of the mixed V-line transforms acting on symmetric \(m\)-tensor fields. Based on a suitable decomposition of tensor fields, we derive new inversion formulas for reconstructing special classes of symmetric \(m\)-tensor fields from the corresponding V-line transform data. These theoretical results extend the current understanding of V-line tensor tomography and provide new insights into the structure of these integral transforms.

In addition, we provide a detailed numerical verification and validation of the inversion algorithms corresponding to the theoretical reconstruction formulas established in [1], [2] for vector field and symmetric \(2\)-tensor fields. The performance of the reconstruction procedures is examined using a variety of representative phantoms, and all numerical experiments are carried out in MATLAB. Furthermore, reconstructions are performed under different levels of additive noise to assess the practical performance of the proposed algorithms.

The remainder of this article is organized as follows. In Section 2, we introduce the straight line and V-line transforms considered in this work together with the notation used throughout the article. Section 3 presents the main theoretical results of this article. Section 4 is devoted to a detailed proof of the main results. Section 5 describes the phantoms used in the numerical experiments and explains the procedure for generating the forward data. This section is further divided into two subsections, Subsections 5.1 and Subsections 5.2, which focus on the reconstruction of vector fields and symmetric \(2\)-tensor fields, respectively. Finally, we conclude the article with acknowledgements in Section 6.

2 Definition, Notations and Preliminaries↩︎

This introductory section presents the notation and definitions used throughout this article. The regular fonts are used to denote scalars or scalar-valued functions (for example, \(x_1\), \(x_2\), \(f\), and \(g\), etc), whereas bold fonts are used to represent vectors, vector fields, or tensor fields in \(\mathbb{R}^2\) (for example, \(\boldsymbol{\textit{f}}\), \(\boldsymbol{\textit{x}}\), \(\boldsymbol{\textit{v}}\), etc).
The disk of radius \(R\) centered at the origin is denoted by \(\mathbb{D}_R\), and its boundary is represented by \(\partial \mathbb{D}_R\). We denote \(C_0^\infty(\mathbb{D}_R, S^m(\mathbb{R}^2))\) to be the space of the symmetric \(m\)-tensor fields whose components belong to \(C_0^\infty(\mathbb{D}_R)\), the space of smooth functions compactly supported in \(\mathbb{D}_R\). Let \(\omega^{1}, \omega^{2},\dots, \omega^{m}\) are vectors in \(\mathbb{R}^2\), then the tensor product \(\omega^{1}\otimes \omega^{2}\otimes \dots\otimes \omega^{m}\) is defined component-wise by \[\begin{align} \left( \omega^{1}\otimes \omega^{2}\otimes \dots\otimes \omega^{m}\right)_{i_1i_2\cdots i_m} = \omega^{1}_{i_1}\omega^{2}_{i_2}\dots \omega^{m}_{i_m}, \qquad i_1,i_2,\cdots,i_m \in \{1, 2\}. \end{align}\] The corresponding symmetric tensor product is defined by \[\begin{align} \omega^{1}\omega^{2}\dots \omega^{m} = \frac{1}{m!}\sum\limits_{\sigma\in \Pi_{m}} \omega^{\sigma(1)}\otimes \omega^{\sigma(2)}\otimes \dots\otimes \omega^{\sigma(m)}, \end{align}\] where \(\Pi_m\) denotes the symmetric group consisting of all permutations of the set \(\{1,2,\ldots,m\}\), for more deetails see [37]. The standard inner product on \(S^{m}(\mathbb{R}^{2})\) is defined by \[\begin{align} \left\langle \boldsymbol{\textit{f}},\boldsymbol{\textit{g}}\right\rangle = f_{i_1i_2\cdots i_m}g^{i_1i_2\cdots i_m} \end{align}\] and the corresponding norm is denoted by \(|\cdot|\). Throughout this article, we use the “Einstein summation convention” according to which summation over repeated indices is understood automatically from \(1\) to \(2\). F or \(\boldsymbol{\textit{u}}=(u^1,u^2)\in\mathbb{R}^2\), the notation \(\boldsymbol{\textit{u}}^m\) represents the symmetric tensor obtained by taking the \(m\)-fold symmetric product of \(\boldsymbol{\textit{u}}\) with itself. Its components are given by \[\begin{align} (\boldsymbol{\textit{u}}^m)_{i_1\cdots i_m} = u^{i_1}\cdots u^{i_m}, \qquad i_1,\ldots,i_m\in\{1,2\}. \end{align}\] Thus, we have \(\boldsymbol{\textit{u}}^m \in S^{m}(\mathbb{R}^2)\).
Next, we introduce the differential operators used later in this article. Compositions of the classical differential operators, namely the gradient and orthogonal gradient operators, naturally give rise to the inner differentiation and inner orthogonal differentiation operators acting on symmetric \(m\)-tensor fields \(\boldsymbol{\textit{f}}\). We also define the divergence and orthogonal divergence operators for symmetric \(m\)-tensor fields. These operators are given as follows:

• The inner differentiation and inner orthogonal differentiation operators \(\mathrm{d},\mathrm{d}^\perp :C_0^\infty(\mathbb{D}_R,S^m(\mathbb{R}^2)) \rightarrow C_0^\infty\left(\mathbb{D}_R, S^{m+1}(\mathbb{R}^2)\right)\) are defined as follows: \[\begin{align} ( \mathrm{d}\boldsymbol{\textit{f}})_{i_1\dots i_mj} &:= \frac{1}{m+1}\left(\frac{\partial{f_{i_1\dots i_m}}}{\partial x_j}+\sum\limits_{k=1}^m\frac{\partial{f_{i_1\dots i_{k-1}ji_{k+1}\dots i_m}}}{\partial x_{i_k}}\right) \quadand\\ (\mathrm{d}^\perp\boldsymbol{\textit{f}})_{i_1\dots i_mj} &:= \frac{1}{m+1}\left((-1)^j\frac{\partial{f_{i_1\dots i_m}}}{\partial x_{3-j}}+\sum\limits_{k=1}^m(-1)^{i_k}\frac{\partial{f_{i_1\dots i_{k-1}ji_{k+1}\dots i_m}}}{\partial x_{3-i_k}}\right). \end{align}\]

• The divergence and the orthogonal divergence operators \(\delta, \delta^\perp:C_0^\infty\left(\mathbb{D}_R,S^{m}(\mathbb{R}^2)\right) \rightarrow C_0^\infty\left(\mathbb{D}_R,S^{m-1}(\mathbb{R}^2)\right)\) are defined as follows: \[\begin{align} (\delta \boldsymbol{\textit{f}})_{i_1\dots i_{m-1}} &:= \frac{\partial{f_{i_1\dots i_{m-1}1}}}{\partial x_1}+\frac{\partial{f_{i_1\dots i_{m-1}2}}}{\partial x_2} \quadand\\(\delta^\perp \boldsymbol{\textit{f}})_{i_1\dots i_{m-1}} &:= \frac{\partial{f_{i_1\dots i_{m-1}2}}}{\partial x_1} - \frac{\partial{f_{i_1\dots i_{m-1}1}}}{\partial x_2}. \end{align}\]

We now recall the definitions of the V-line (or broken ray) and the corresponding V-line transforms for symmetric \(m\)-tensor fields, which were previously studied in [1], [2], [33]. For the definitions and background on straight-line transforms, we refer the reader to [37], [38].

More precisely, \(BR(\beta, d)\) is defined by: \[\begin{align} \label{eq:definition32of32BR40beta44d41} BR(\beta,d) = \left\{\boldsymbol{\textit{x}}_\beta + t \boldsymbol{\textit{u}}_\beta: 0\leq t \leq d\right\} \cup \left\{\boldsymbol{\textit{x}}_\beta + d \boldsymbol{\textit{u}}_\beta + s \boldsymbol{\textit{v}}_\beta: 0 \leq s < \infty\right\}. \end{align}\tag{1}\]

Here \(\boldsymbol{\textit{x}}_\beta = (R \cos \beta, R \sin \beta)\) denotes the point on the boundary \(\partial\mathbb{D}_R\) from which the broken ray originates. The vectors \(\boldsymbol{\textit{u}}_\beta = -(\cos \beta, \sin \beta)\) and \(\boldsymbol{\textit{v}}_\beta = -\left(\cos (\theta + \beta), \sin (\theta +\beta)\right)\) are unit vectors representing the directions of the first and second segments of the V-line, respectively. Their corresponding orthogonal unit vectors are given by \(\boldsymbol{\textit{u}}_{\beta}^\perp = (\sin \beta, - \cos \beta)\) and \(\boldsymbol{\textit{v}}_\beta^{\perp} = \left(\sin (\theta +\beta), -\cos (\theta + \beta)\right)\) respectively. For each pair \((\beta,d)\), the unit vectors \(\boldsymbol{\textit{u}}_\beta\) and \(\boldsymbol{\textit{v}}_\beta\) are uniquely determined and will be used throughout the subsequent definitions of the V-line transforms acting on symmetric \(m\)-tensor fields.

Subsequently, we also introduce the straight-line ray transforms for the symmetric \(m\)-tensor field, which will be used multiple times later in this article. These transforms are defined by integrating tensor fields along straight lines in \(\mathbb{R}^{2}\). For \(\psi\in[0,2\pi)\) and \(p\in\mathbb{R}\), let \[\begin{align} L(\psi,p):=\{(y_1,y_2)\in\mathbb{R}^{2}:y_1\cos\psi+y_2\sin\psi=p\} \end{align}\] denote the straight line whose signed distance from the origin is \(p\) and whose normal vector is \(\boldsymbol{\textit{w}}=(\cos\psi,\sin\psi)\). The corresponding unit tangent vector is \(\boldsymbol{\textit{w}}^\perp=(-\sin\psi,\cos\psi),\) which is orthogonal to \(\boldsymbol{\textit{w}}\). The line \(L(\psi,p)\) can therefore be parameterized as \(\{\,p\boldsymbol{\textit{w}}+s\boldsymbol{\textit{w}}^\perp:\; s\in\mathbb{R}\,\}\).

Finally, we state the Radon transform and its inversion.

The injectivity of the Radon transform implies that \(f\) can be reconstructed uniquely from its Radon data. An explicit inversion formula is given by (see [39]): \[\begin{align} \label{eq:Radon39s32inversion32formula} f(\boldsymbol{\textit{x}}) = \frac{1}{2\pi} \left(-\Delta\right)^{1/2} \int_0^{2\pi}\mathcal{R}f((\cos \alpha, \sin \alpha), x \cos \alpha + y \sin \alpha) d\alpha. \end{align}\tag{5}\] Next, we present the decomposition results for symmetric \(m\)-tensor fields that will be needed later.

3 Theoretical Results↩︎

This section is devoted to the statement of the main results of the article. We start by presenting the kernel characterisation of the mixed V-line transforms for symmetric \(m\)-tensor fields. We show that each of these transforms possesses a non-trivial kernel. This behavior is also analogous to the case of the straight-line integral transforms of symmetric \(m\)-tensor fields, where the corresponding kernel characterizations are known to be exact. The next theorem describes the kernel structure of the V-line transforms associated with symmetric \(m\)-tensor fields.

Next, we recall a reconstruction result for symmetric \(m\)-tensor fields from V-line transform data, previously obtained in [1], [2]. The numerical implementation, verification, and validation of the associated inversion algorithms for vector fields and symmetric \(2\)-tensor fields are presented in Section 5.

The next result provides the inversion formulas for the V-line transforms associated with special classes of symmetric \(m\)-tensor fields. More precisely, we consider tensor fields of the form \[\begin{align} \boldsymbol{\textit{f}}=(\mathrm{d}^\perp)^{m-k}\mathrm{d}^k\Psi^{(k)}, \qquad 0\leq k\leq m, \end{align}\] where \(\Psi^{(k)}\) is a scalar potential. For this particular class of symmetric \(m\)-tensor fields, complete reconstruction can be achieved using only the \(k\)-th mixed V-line transform data. This is due to the fact that all remaining components of the tensor field belong to the kernel of the corresponding \(k\)-th mixed V-line transform, which was established in Theorem 2. Consequently, for this special class of symmetric tensor fields, a single mixed V-line transform contains sufficient information for complete recovery.

To prove this result, our approach converts the information obtained from the V-line transform data into corresponding straight-line transform data. Once this connection is established, the inversion formulas for straight-line transforms can be applied to reconstruct the unknown vector field or symmetric \(m\)-tensor field.

It is important to note that all these theorems require certain support restrictions on \(\boldsymbol{\textit{f}}\), which depend on the fixed scattering angle \(\theta\). These restrictions arise naturally from the reconstruction method employed in the proofs. In particular, as illustrated in Figure 1 (b), straight-line transform data corresponding to lines outside the disk of radius \(R\sin{\theta}\) cannot be generated from the available V-line data. This reconstruction strategy was originally introduced by Gaik in [33] for scalar functions.

The detailed proofs of Theorem 3 for the vector field and the symmetric \(m\)-tensor field can be found in [1], [2]. The results stated in Theorem 2 and Theorem 4 are new. Detailed proofs of both theorems are provided in the Section 4.

4 Proofs of the Main Results↩︎

This section is devoted to the proofs of Theorems 3.1 and 3.3. These results establish the kernel characterization of the mixed V-line transforms acting on symmetric \(m\)-tensor fields and derive inversion formulas for the reconstruction of special classes of symmetric \(m\)-tensor fields from the corresponding V-line transform data.

4.1 Proof of Theorem 2↩︎

4.2 Proof of Theorem 4↩︎

5 Numerical Implementation↩︎

This section is devoted entirely to the numerical implementation of the proposed method and to evaluating its performance on different phantoms. We present the numerical reconstruction only for vector fields \((m=1)\) and symmetric 2-tensor fields \((m=2)\). To avoid repetition, we have used different combinations of data for different phantoms and vector/tensor fields. Because of this, our numerical validation is primarily divided into two parts:

1. First, we generate the forward data for a given field \(\boldsymbol{\textit{f}}\) over the domain of \([-1, 1]\times [-1, 1]\). For the vector field case, we generate V-line data, longitudinal V-line transform \(\mathcal{L}\boldsymbol{\textit{f}}\) and transverse V-line transform \(\mathcal{T}\boldsymbol{\textit{f}}\), and for the tensor field case, we generate V-line data, longitudinal V-line transform \(\mathcal{L}\boldsymbol{\textit{f}}\), transverse V-line transform \(\mathcal{T}\boldsymbol{\textit{f}}\) and mixed V-line transform \(\mathcal{M}\boldsymbol{\textit{f}}\).

2. Utilizing the data generated in the first step, we reconstruct the underlying components of the vector/tensor fields. More specifically, we numerically reconstruct the vector field \(\boldsymbol{\textit{f}}= (f_1, f_2)\) using the data obtained from the V-line data, \(\mathcal{L}\boldsymbol{\textit{f}}\) and \(\mathcal{T}\boldsymbol{\textit{f}}\). Analogously, for symmetric \(2\)-tensor field \(\boldsymbol{\textit{f}}= \begin{pmatrix} f_{11} & f_{12}\\ f_{21} & f_{22} \end{pmatrix}\), we use the information of V-line data \(\mathcal{L}\boldsymbol{\textit{f}}\), \(\mathcal{T}\boldsymbol{\textit{f}}\), and \(\mathcal{M}\boldsymbol{\textit{f}}\). In this case, we will identify \(\boldsymbol{\textit{f}}\) with \((f_{11}, f_{12}, f_{22})\), since \(\boldsymbol{\textit{f}}\) is symmetric.

A more detailed discussion of both the data generation process and the subsequent numerical reconstruction begins with the vector field case in Subsection 5.1, and then the tensor field case in Subsection 5.2.

5.1 Reconstruction of Vector Fields↩︎

5.1.1 Data Generation↩︎

To generate data for the numerical simulations, we consider a circular domain of radius \(R=1\) and fix the scattering angle at \(\theta = \pi/4\). To satisfy the theoretical requirement that the field be supported within a disk of radius \(R \sin \theta\), a spatial mask is applied to the considered phantom. We compute the discrete forward data by evaluating the longitudinal (\(\mathcal{L}\boldsymbol{\textit{f}}\)) and transverse (\(\mathcal{T}\boldsymbol{\textit{f}}\)) V-line transforms over a grid of source angles \(\beta\) and scattering depths \(d\). For each pair \((\beta, d)\), the numerical integration is performed in two stages:

1. Incident Segment: The ray is traced from the boundary source to the vertex at depth \(d\) along the direction \(\boldsymbol{\textit{u}}\). Since discrete points along this path rarely align with the Cartesian grid, 2D bilinear interpolation is used to sample the vector field components \((f_1, f_2)\). We then compute the inner products of \((f_1, f_2)\) with \(\boldsymbol{\textit{u}}\) and \(\boldsymbol{\textit{u}}^{\perp}\).

2. Scattered Segment: The path continues from the vertex along the deflected direction \(\boldsymbol{\textit{v}}\). We follow the same interpolation and integration procedure to compute the contribution coming from this branch.

The line integrals are approximated using the composite trapezoidal rule. Finally, the total discrete measurements for \(\mathcal{L}\boldsymbol{\textit{f}}(\beta,d)\) and \(\mathcal{T}\boldsymbol{\textit{f}}(\beta,d)\) are obtained by summing the contributions from both segments.

5.1.2 Numerical Reconstruction↩︎

In Figures 2–9, we demonstrate the recovery of the scalar functions \(f_1\) and \(f_2\), components of a vector field \(\boldsymbol{\textit{f}}\) across various phantoms, both with noise and without noise. To achieve this, we first utilize our V-line transform data \(\mathcal{L}\boldsymbol{\textit{f}}\) and \(\mathcal{T}\boldsymbol{\textit{f}}\) to generate longitudinal and transverse ray transforms, \(\mathcal{I}\boldsymbol{\textit{f}}\) and \(\mathcal{J}\boldsymbol{\textit{f}}\) (straight line transforms), via the following relation (for more details see [1]): \[\begin{align} \label{relation32btw32V-line32and32st32line} \begin{cases} \mathcal{I}\boldsymbol{\textit{f}}(\psi_{\beta},t_{d}) = \mathcal{L}\boldsymbol{\textit{f}}(\beta,d) - \mathcal{L}\boldsymbol{\textit{f}}(\beta + \pi,2R-d) - \mathcal{L}\boldsymbol{\textit{f}}(\beta, 2R),\\[5pt] \mathcal{J}\boldsymbol{\textit{f}}(\psi_{\beta},t_{d}) = \mathcal{T}\boldsymbol{\textit{f}}(\beta,d) - \mathcal{T}\boldsymbol{\textit{f}}(\beta + \pi,2R-d) - \mathcal{T}\boldsymbol{\textit{f}}(\beta, 2R), \end{cases} \end{align}\tag{15}\] where \(d \in [0,2R]\), \(\beta \in [0, 2\pi)\), \(\psi_{\beta} = \beta + \theta + \pi/2\), and \(t_{d} = (R-d)\sin(\pi + \theta)\).
Recall, for \(\psi \in [0, 2\pi), \;p \in \mathbb{R},\) the longitudinal ray transform \(\mathcal{I}\boldsymbol{\textit{f}}\) and transverse ray transform \(\mathcal{J}\boldsymbol{\textit{f}}\) are given by: \[\label{eq:def32LRT} \mathcal{I}\boldsymbol{\textit{f}}(\psi,p) = \mathcal{I}\boldsymbol{\textit{f}}(\boldsymbol{\textit{w}},p) := \int_\mathbb{R} \boldsymbol{\textit{w}}^\perp \cdot \boldsymbol{\textit{f}}(p \boldsymbol{\textit{w}}+ s\boldsymbol{\textit{w}}^\perp)\,ds,\tag{16}\] \[\label{eq:def32TRT} \mathcal{J}\boldsymbol{\textit{f}}(\psi,p) = \mathcal{J}\boldsymbol{\textit{f}}(\boldsymbol{\textit{w}},p) := \int_\mathbb{R} \boldsymbol{\textit{w}}\cdot \boldsymbol{\textit{f}}(p \boldsymbol{\textit{w}}+ s\boldsymbol{\textit{w}}^\perp)\,ds.\tag{17}\] These transforms can be cast in matrix form to isolate the Radon transforms of the individual scalar components: \[\label{eq:matrix95system} \begin{bmatrix} \mathcal{I}\boldsymbol{\textit{f}}\\ \mathcal{J}\boldsymbol{\textit{f}} \end{bmatrix}(\boldsymbol{\textit{w}}, p) = \begin{bmatrix} -w_2 & w_1 \\ w_1 & w_2 \end{bmatrix} \begin{bmatrix} \mathcal{R}f_1 \\ \mathcal{R}f_2 \end{bmatrix}(\boldsymbol{\textit{w}}, p), \quadhere\boldsymbol{\textit{w}}= (w_1, w_2) = (\cos \psi, \sin \psi).\tag{18}\] Because the transformation matrix in 18 is invertible, we can point-wise solve this system to isolate the component-wise Radon transforms, \(\mathcal{R}f_1\) and \(\mathcal{R}f_2\). Finally, applying the inverse Radon transform yields the reconstructed scalar components \(f_1\) and \(f_2\).

We present numerical experiments starting with a baseline, noise-free case. We evaluate the performance of our inversion algorithm using three distinct phantom classes defined on the domain \([-1, 1] \times [-1, 1]\). Each class introduces varying degrees of mathematical regularity and geometric complexity: smooth bump functions (Phantom 1 - PH1), characteristic functions of overlapping disks (Phantom 2 - PH2), and characteristic functions of rectangular annuli (Phantom 3 - PH3). As a first test case, we define the components \(f_1\) and \(f_2\) using smooth Gaussian bump functions centered at different locations with distinct supports. The mathematical function for each bump is defined as: \[f_i(x,y) = \begin{cases} e^{-s^2 / \left(s^2 - \left[(x - a)^2+(y - b)^2\right]\right)}, & (x - a)^2+(y - b)^2 < s^2, \\[2pt] 0, & (x - a)^2+(y - b)^2 \geq s^2, \end{cases}\] where \((a, b)\) denotes the center of the bump function, and \(s\) represents the radius of its support. The selected parameters are summarized in Table 1. We will refer to this phantom as PH1. Later, for the tensor field, there will be one more component, which we will also call PH1.

Figure 2 displays the results for this smooth test case. The first column shows the ground-truth scalar components \(f_1\) and \(f_2\). The second column illustrates the simulated forward data \(\mathcal{L}\boldsymbol{\textit{f}}(\beta, d)\) and \(\mathcal{T}\boldsymbol{\textit{f}}(\beta,d)\) plotted against the angular variable \(\beta\) (in degrees) and the distance parameter \(d\). The final column displays the reconstructed components obtained via the proposed inversion method. We see that the reconstructed scalar components are really good in the absence of noise (see Table 4 for the relative error).

To investigate the stability of the proposed inversion method, we next introduce varying levels of Gaussian noise to the forward projection data. Figure 3 illustrates the reconstruction of phantom 1 under noise levels of \(5\%, 10\%\) and \(20\%\). In the reconstruction, we observe that increasing the noise level introduces visible granular artifacts into the spatial domain. However, we observe that the locations, shapes, and amplitudes of \(f_1\) and \(f_2\) are well reconstructed. Thus, the reconstruction remains robust in the presence of noise.

We now consider another phantom, which we refer to as Ph2. In this class, each component of the phantom is represented as a weighted combination of three characteristic functions of disks having different radii \(r\) and centers \((a,b)\). The values of these parameters are given in Table 2 :

From Figure 4, it is evident that the method handles the sharp discontinuities of the phantom with minimal visual artifacts.

Extending the validation of the inversion scheme, Figure 5 evaluates the numerical stability of the algorithm against data perturbation by reconstructing the piecewise constant field from noisy projection. As in the previous cases, the leftmost column shows the original components (\(f_1\) and \(f_2\)), while the subsequent columns show the reconstructed fields after corrupting the transform data with \(5\%, 10\%,\) and \(20\%\) noise, respectively. We observe that reconstructed components have prominent granular variations and distinct circular boundary artifacts. But the numerical implementation also shows robustness, as it can capture the fundamental features even under severe \(20 \%\) noise conditions.

Further expanding on the validation of the inversion framework, we consider a vector field characterized by non-convex, piecewise-constant geometries, specifically rotated square-ring structures. This last class of phantoms consists of square annuli with identical dimensions, centered at different locations, and rotated at different angles with respect to the \(x\)-axis. We will refer to this phantom as Ph3. For details about the locations and rotations, please refer to the Table 3.

Table 4 summarizes the relative errors in reconstructing \(f_1\) and \(f_2\) from V-line data \(\mathcal{L} \boldsymbol{\textit{f}}\) and \(\mathcal{T}\boldsymbol{\textit{f}}\) for the three phantoms (PH1–PH3) under different noise levels. The results show that the reconstruction error increases as the noise level rises. Figure 6 demonstrates the algorithm’s ability to accurately reconstruct the components of a vector field defined on the phantom, which have sharp corners and hollow interiors.

Further, Figure 7 shows the successful reconstruction of the non-convex square ring phantom. This figure evaluates the algorithm’s stability by introducing \(5 \%, 10 \%\) and \(20 \%\) noise to the projected data.

Table 4 summarizes the relative errors in reconstructing \(f_1\) and \(f_2\) from V-line data \(\mathcal{L} \boldsymbol{\textit{f}}\) and \(\mathcal{T}\boldsymbol{\textit{f}}\) for the three phantoms (PH1–PH3) under different noise levels. The results show that the reconstruction error increases as the noise level rises. Next, we present a brief comparison of reconstructions from different breaking angles \(\theta\). Recall that this angle \(\theta\) directly dictates the observable field of view (FOV), as seen from the theoretical results (Theorem 3) that the reconstruction works if the unknown field is supported in a disk of radius \(R \sin(\theta)\). Consequently, small angles (e.g., \(15^\circ\)) yield a severely restricted FOV; any part of the vector field outside this small central region suffers from severe truncation artifacts. In contrast, as the opening angle widens (e.g., \(75^\circ\)), the valid mathematical FOV expands to encompass nearly the entire domain, allowing full reconstruction. We have presented this comparison for Phantom 2 (see Figure 8) and Phantom 3 (see Figure 9) with \(\theta\) taking the values \(15^\circ\), \(30^\circ\), \(45^\circ\), \(60^\circ,\) and \(75^\circ\).

5.2 Reconstruction of Tensor Fields↩︎

To generate the data, we will follow the same procedure as discussed in subsection 5.1. The only difference is that we will now consider the three components of \(\boldsymbol{\textit{f}}=(f_{11}, f_{12},f_{22})\), and we compute \(\mathcal{L}\boldsymbol{\textit{f}}\), \(\mathcal{M}\boldsymbol{\textit{f}}\), and \(\mathcal{T}\boldsymbol{\textit{f}}\). Now the goal is to recover \(f_{11}, f_{12}\) and \(f_{22}\) from the given \(\mathcal{L}\boldsymbol{\textit{f}}\), \(\mathcal{M}\boldsymbol{\textit{f}}\), and \(\mathcal{T}\boldsymbol{\textit{f}}\). Figures [10-13], demonstrates the recovery of the components \(f_{11}\), \(f_{12}\), and \(f_{22}\) of a tensor field \(\boldsymbol{\textit{f}}\) across various phantoms, both with and without noise. To achieve this, we follow a similar approach analogous to the vector field case, i.e., we first compute the following formulas to get \(\mathcal{I}\boldsymbol{\textit{f}}\), \(\mathcal{K}\boldsymbol{\textit{f}}\), and \(\mathcal{J}\boldsymbol{\textit{f}}\) using \(\mathcal{L}\boldsymbol{\textit{f}}\), \(\mathcal{M}\boldsymbol{\textit{f}}\), and \(\mathcal{T}\boldsymbol{\textit{f}}\): \[\begin{align} \label{eq59:relation32btw32V-line32and32st32line32transform32for32tensor32field} \begin{cases} \mathcal{I}\boldsymbol{\textit{f}}{{(\psi_{\beta},t_{d})}} = \mathcal{L}\boldsymbol{\textit{f}}(\beta,d) + \mathcal{L}\boldsymbol{\textit{f}}(\beta + \pi,2R-d) - \mathcal{L}\boldsymbol{\textit{f}}(\beta, 2R),\\[4pt] \mathcal{K}\boldsymbol{\textit{f}}(\psi_{\beta},t_{d}) = \mathcal{M}\boldsymbol{\textit{f}}(\beta,d) -\mathcal{M}\boldsymbol{\textit{f}}(\beta + \pi,2R-d) - \mathcal{M}\boldsymbol{\textit{f}}({\beta,2R}),\\[4pt] \mathcal{J}\boldsymbol{\textit{f}}(\psi_{\beta},t_{d}) = \mathcal{T}\boldsymbol{\textit{f}}(\beta,d) + \mathcal{T}\boldsymbol{\textit{f}}(\beta + \pi,2R-d) - \mathcal{T}\boldsymbol{\textit{f}}(\beta, 2R), \end{cases} \end{align}\tag{19}\] where \(d\in[0,2R]\), \(\beta \in [0, 2\pi)\), \(\psi_{\beta} = \beta + \theta + \pi/2 \; { and } \; t_{d} = (R-d)\sin(\pi + \theta)\). For \(\psi \in [0, 2\pi) { and } p \in \mathbb{R}\), we recall the longitudinal ray transform \(\mathcal{I}\boldsymbol{\textit{f}}\), mixed ray transform \(\mathcal{K}\boldsymbol{\textit{f}}\) and transverse ray transform \(\mathcal{J}\boldsymbol{\textit{f}}\), which are defined by \[\label{eq:def32LRT32for32tensor} \mathcal{I}\boldsymbol{\textit{f}}(\psi,p) = \mathcal{I}\boldsymbol{\textit{f}}(\boldsymbol{\textit{w}},p) := \int_\mathbb{R} \left\langle\boldsymbol{\textit{f}}(p \boldsymbol{\textit{w}}+ s\boldsymbol{\textit{w}}^\perp), (\boldsymbol{\textit{w}}^\perp)^{2} \right\rangle\,ds,\tag{20}\] \[\label{eq:def32MRT} \mathcal{K}\boldsymbol{\textit{f}}(\psi,p) = \mathcal{K}\boldsymbol{\textit{f}}(\boldsymbol{\textit{w}},p) := \int_\mathbb{R} \left\langle \boldsymbol{\textit{f}}(p \boldsymbol{\textit{w}}+ s\boldsymbol{\textit{w}}^\perp), \boldsymbol{\textit{w}}^{\perp}\odot\boldsymbol{\textit{w}}\right\rangle\,ds,\tag{21}\] \[\label{eq:def32TVT32for32tensor} \mathcal{J}\boldsymbol{\textit{f}}(\psi,p) = \mathcal{J}\boldsymbol{\textit{f}}(\boldsymbol{\textit{w}},p) := \int_\mathbb{R} \left\langle \boldsymbol{\textit{f}}(p \boldsymbol{\textit{w}}+ s\boldsymbol{\textit{w}}^\perp), \boldsymbol{\textit{w}}^{2} \right\rangle\,ds,\tag{22}\] respectively. Observe that, \[\begin{align} \begin{bmatrix} \mathcal{I} \boldsymbol{\textit{f}}\\ \mathcal{K}\boldsymbol{\textit{f}}\\ \mathcal{J}\boldsymbol{\textit{f}} \end{bmatrix} = \overbrace{ \begin{bmatrix} w_2^2 & - \frac{w_1 w_2 + w_2 w_1}{2} & w_1^2 \\ -w_2 w_1 & \frac{w_!^2 - w_2 ^2}{2} & w_! w_2 \\ w_1^2 & \frac{w_1 w_2 + w_1 w_2}{2} & w_2^2 \end{bmatrix}}^{W} \begin{bmatrix} \mathcal{R}f_{11} \\ \mathcal{R}f_{12} \\ \mathcal{R}f_{22} \end{bmatrix} \end{align}\] Note that matrix \(W\) is invertible, so we can solve this system of equations to get \(\mathcal{R} f_{11}\), \(\mathcal{R} f_{12}\), and \(\mathcal{R} f_{22}\), and further applying the inverse Radon transform, we can get the components of \(\boldsymbol{\textit{f}}\).

To numerically validate the results, we implement the algorithm as in the earlier vector field case. To avoid repetition, we show the reconstruction for the phantom Ph2 and Ph3, with and without noise. Other cases can also be handled analogously. The parameters for these phantoms are given in Table 5 and 6.

Figure 10 and 12 shows the reconstruction of the \(f_{11}, f_{12}\) and \(f_{22}\) for Ph2 and Ph3 in the absence of noise. This shows that the algorithm can reconstruct sharp discontinuities for both convex and non-convex geometries with almost no artifacts. Figure 11 and 13 show that we obtain a nice reconstruction of the components of \(\boldsymbol{\textit{f}}\) in the presence of \(5 \%, 10 \%\) and \(20 \%\) noise. There are sufficient granular artifacts but it can capture the basic features despite having \(20 \%\) noise in the data. Table 7 gives the relative error for various cases, including PH2 and PH3 with and without noise.

6 Acknowledgements↩︎

Rahul Bhardwaj gratefully acknowledges the partial financial support provided by the FIST programme of the Department of Science and Technology (DST), Government of India, under Grant No. SR/FST/MS-I/2018/22(C).

Data availability statement. No datasets were generated or analyzed during the current study; therefore, data sharing is not applicable.

Conflict of interest. The authors declare that they have no conflicts of interest regarding the research, authorship, and/or publication of this article.

References↩︎

---

1. Department of Mathematics, Indian Institute of Technology, Ropar, Punjab - 140001, India. bhardwaj161067@gmail.com↩︎

2. Department of Mathematics, Indian Institute of Technology, Gandhinagar, Gujarat - 382355, India. madhu.gupta@iitgn.ac.in↩︎