Torus computed tomography

TL;DR

This paper introduces Torus CT, a novel 2D computed tomography method based on the geometry of the flat torus, providing new inversion formulas, regularization strategies, and promising simulation results.

Contribution

The paper develops a new CT inversion technique using torus geometry, including formulas, regularization, and implementation details, advancing the field of tomographic reconstruction.

Findings

New inversion formulas for the Radon transform on the flat torus

Regularization strategy with stability estimates

Successful implementation with simulated data

Abstract

We present a new computed tomography (CT) method for inverting the Radon transform in 2D. The idea relies on the geometry of the flat torus, hence we call the new method Torus CT. We prove new inversion formulas for integrable functions, solve a minimization problem associated to Tikhonov regularization in Sobolev spaces and prove that the solution operator provides an admissible regularization strategy with a quantitative stability estimate. This regularization is a simple post-processing low-pass filter for the Fourier series of a phantom. We also study the adjoint and the normal operator of the X-ray transform on the flat torus. The X-ray transform is unitary on the flat torus. We have implemented the Torus CT method using Matlab and tested it with simulated data with promising results. The inversion method is meshless in the sense that it gives out a closed form function that can be evaluated at any point of interest.

Tables

Table 1. Table 1 . The regularization parameters ( α , s ) 𝛼 𝑠 (\alpha,s) that give the best reconstructions with respect to the L p superscript 𝐿 𝑝 L^{p} -norms with p = 1 , 2 , ∞ 𝑝 1 2 p=1,2,\infty ; respective error  ϵ p α , s superscript subscript italic-ϵ 𝑝 𝛼 𝑠 \epsilon_{p}^{\alpha,s} ; and error  ϵ p 0 , 0 superscript subscript italic-ϵ 𝑝 0 0 \epsilon_{p}^{0,0} of non-regularized reconstruction.

norm	Shepp–Logan				Flag
$p$	$\alpha$	$s$	${ϵ}_p^{\alpha ,s}$	${ϵ}_p^{0,0}$	$\alpha$	$s$	${ϵ}_p^{\alpha ,s}$	${ϵ}_p^{0,0}$
1	$0.050$	$0.69$	$62$%	$112$%	$0.025$	$0.71$	$41$%	$69$%
2	$0.025$	$0.61$	$48$%	$70$%	$0.025$	$0.68$	$29$%	$45$%
$\infty$	$0.025$	$0.56$	$75$%	$112$%	$0.025$	$0.78$	$73$%	$106$%

Table 2. Table 2 . Errors in reconstructions computed with 𝒜 𝕋 2 subscript 𝒜 superscript 𝕋 2 \mathcal{A}_{\mathbb{T}^{2}} and the FBP.

	Shepp–Logan	Flag	Rotated Flag	Shepp–Logan	Flag	Rotated Flag
	${𝒜}_{{𝕋}^2}$ with noiseless data ($ℰ=0$)			${𝒜}_{{𝕋}^2}$ with noisy data ($ℰ\sim 𝒩(0,{\sigma }^2)$)
${ϵ}_1^{0,0}$	$313$%	$298$%	$302$%	$310$%	$300$%	$301$%
${ϵ}_2^{0,0}$	$161$%	$171$%	$168$%	$162$%	$172$%	$167$%
${ϵ}_{\infty }^{0,0}$	$75$%	$108$%	$121$%	$79$%	$125$%	$125$%
				Regularized reconstruction from noisy data
${ϵ}_1^{\alpha ,s}$				$331$%	$305$%	$303$%
${ϵ}_2^{\alpha ,s}$				$170$%	$174$%	$173$%
${ϵ}_{\infty }^{\alpha ,s}$				$53$%	$99$%	$100$%
	FBP with torus optimized angles			FBP with evenly distributed angles
${ϵ}_1$	$73$%	$59$%	$67$%	$64$%	$55$%	$56$%
${ϵ}_2$	$59$%	$41$%	$51$%	$54$%	$45$%	$45$%
${ϵ}_{\infty }$	$155$%	$87$%	$127$%	$129$%	$93$%	$120$%

Table 3. Table 3 . Reconstruction errors  ϵ p 0 , 0 superscript subscript italic-ϵ 𝑝 0 0 \epsilon_{p}^{0,0} of (non-regularized) reconstructions from rotational data sets.

Shepp–Logan	${\Theta }_1^{\mathrm{SL}}$	${\Theta }_2^{\mathrm{SL}}$	${\Theta }_3^{\mathrm{SL}}$	${\Theta }_4^{\mathrm{SL}}$	${\Theta }_5^{\mathrm{SL}}$	${\Theta }_6^{\mathrm{SL}}$	${\Theta }_7^{\mathrm{SL}}$	${\Theta }_8^{\mathrm{SL}}$	${\Theta }_9^{\mathrm{SL}}$
$p=1$	$113$%	$82$%	$69$%	$61$%	$56$%	$52$%	$49$%	$46$%	$44$%
$p=2$	$72$%	$54$%	$46$%	$42$%	$40$%	$38$%	$36$%	$34$%	$33$%
$p=\infty$	$97$%	$89$%	$84$%	$77$%	$77$%	$73$%	$70$%	$69$%	$68$%
Flag	${\Theta }_1^F$	${\Theta }_2^F$	${\Theta }_3^F$	${\Theta }_4^F$	${\Theta }_5^F$	${\Theta }_6^F$
$p=1$	$22$%	$18$%	$16$%	$14$%	$12$%	$12$%
$p=2$	$25$%	$21$%	$18$%	$16$%	$15$%	$14$%
$p=\infty$	$96$%	$86$%	$79$%	$66$%	$68$%	$69$%

Equations

$$\mathcal{I}f(x,v)={\mathcal{I}}_{v}f(x)={\int }_0^{1}f(x+tv)\mathrm{d}t,x\in {\mathbb{T}}^2,v\in {\mathbb{Z}}^2,$$

$$\hat{f}(k)=\{\begin{array}{l}{\int }_0^1{\mathcal{I}}_{v}f(0,y)\exp (-2\pi i{k}_{2}y)dy,k_2\ne 0\\ {\int }_0^1{\mathcal{I}}_{v}f(x,0)\exp (-2\pi i{k}_{1}x)dx,k_1\ne 0.\end{array}$$

$$\hat{f}(k)={\int }_0^1{\mathcal{I}}_{(1,0)}f(0,y)dy={\int }_0^1{\mathcal{I}}_{(0,1)}f(x,0)dx.$$

$$\underset{f\in H^r({\mathbb{T}}^2)}{\mathrm{argmin}}\left({\Vert \mathcal{I}f-g\Vert }_{H^r({\mathbb{T}}^2\times Q)}^2+\alpha {\Vert f\Vert }_{H^s({\mathbb{T}}^2)}^2\right).$$

$$\underset{\epsilon \to 0}{\lim }\underset{{\Vert g\Vert }_{H^t({\mathbb{T}}^2\times Q)}\le \epsilon }{\sup }{\Vert P_{\alpha (\epsilon )}^s{\mathcal{I}}^{\ast }(\mathcal{I}f+g)-f\Vert }_{H^r({\mathbb{T}}^2)}=0,$$

$${\Vert P_{\alpha }^s{\mathcal{I}}^{\ast }(\mathcal{I}f+g)-f\Vert }_{H^r({\mathbb{T}}^2)}\le {\alpha }^{\delta /2s}C(\delta /2s){\Vert f\Vert }_{H^{r+\delta }({\mathbb{T}}^2)}+\frac{\epsilon }{\alpha },$$

$$\mathcal{I}f(x,v):={\int }_0^{1}f({\gamma }_{x,v}(t))\mathrm{d}t,f\in \mathcal{T},x\in [0,1{]}^{2}v\in {\mathbb{Z}}^2\setminus 0.$$

$$[\mathcal{I}f(\cdot ,v)](\eta ):=(f,\mathcal{I}\eta (\cdot ,v)),\eta \in \mathcal{T}$$

$$f(x)=\sum _{k\in {\mathbb{Z}}^2}\hat{f}(k)e^{2\pi ik\cdot x}$$

$$\widehat{\mathcal{I}f}(k,v)=\{\begin{array}{ll}\hat{f}(k)& k\cdot v=0\\ 0& k\cdot v\ne 0.\end{array}$$

$$(f,\phi )={\int }_{{\mathbb{T}}^2}\tilde{f}\phi dm,\forall \phi \in C^{\infty }({\mathbb{T}}^2).$$

$${\phi }_{v,m}(a,b)=a\frac{v}{\left|v_2\right|}+(0,b)+w_m,a\in [0,\left|\frac{v_2}{v_1}\right|),b\in [0,1).$$

$${\mathcal{I}}_{v}f(p):={\int }_0^{1}f(p+tv)dt\text{for a.e. p∈T2}.$$

$$\hat{f}(k)={\int }_0^1{\int }_0^{1}f(x,y)\exp (-2\pi ik\cdot (x,y))dxdy.$$

$$\hat{f}(k)={\int }_0^1{\int }_0^{|v_2/v_1|}f({\phi }_{v,m}(a,b))\exp \left(-2\pi ik\cdot \left(\frac{a}{\left|v_2\right|}v+(0,b)+w_m\right)\right)\left|\frac{v_1}{v_2}\right|dadb.$$

$$\hat{f}(k)=\left|\frac{v_1}{v_2}\right|{\int }_0^1{\int }_0^{|v_2/v_1|}f\left(a\frac{v}{\left|v_2\right|}+(0,b)+w_m\right)da\exp (-2\pi i{k}_{2}b)db.$$

$$\left|v_1\right|\hat{f}(k)=\left|v_1\right|{\int }_0^1{\mathcal{I}}_{v}f(0,y)\exp (-2\pi i{k}_{2}y)dy.$$

$$\widehat{{\mathcal{I}}_{v}f}(k)=\left|\frac{v_1}{v_2}\right|{\int }_0^1{\int }_0^{|v_2/v_1|}{\mathcal{I}}_{v}f({\phi }_{v,0}(a,b))da\exp (-2\pi i{k}_{2}b)db.$$

$${\mathcal{R}}_{v}f(p)={\int }_{\mathbb{R}}f(p+tv)dt.$$

$${\Vert f\Vert }_{H^s({\mathbb{T}}^2)}^2=\sum _{k\in {\mathbb{Z}}^2}{⟨k⟩}^{2s}{\left|\hat{f}(k)\right|}^2,$$

$${\Vert g\Vert }_{H^s({\mathbb{T}}^2\times Q)}^2={\left|\hat{g}(0,0)\right|}^2+\sum _{k\in {\mathbb{Z}}^2\setminus 0}\sum _{v\in Q}{⟨k⟩}^{2s}{\left|\hat{g}(k,v)\right|}^2$$

$$\widehat{{\mathcal{I}}^{\ast }g}(k)=\hat{g}(k,{\hat{k}}^{\perp }).$$

$$\begin{array}{rl}\left(\mathcal{I}f,g\right)& =\widehat{\mathcal{I}f}(0,0{)}^{\ast }\hat{g}(0,0)+\sum _{k\in {\mathbb{Z}}^2\setminus 0}\sum _{v\in Q}{⟨k⟩}^{2s}\widehat{\mathcal{I}f}(k,v{)}^{\ast }\hat{g}(k,v)\\ & =\sum _{k\in {\mathbb{Z}}^2}{⟨k⟩}^{2s}\hat{f}(k{)}^{\ast }\hat{g}(k,{\hat{k}}^{\perp })\\ & =\left(f,{\mathcal{I}}^{\ast }g\right).\end{array}$$

$${\Vert f\Vert }_{H^s({\mathbb{T}}^2)}^2={\left|\hat{f}(0)\right|}^2+\sum _{v\in Q}\sum _{p\in \mathbb{Z}\setminus 0}{⟨pv⟩}^{2s}{\left|\hat{f}(p{v}^{\perp })\right|}^2$$

$${\Vert \mathcal{I}f-g\Vert }_{H^r({\mathbb{T}}^2\times Q)}^2={\left|\widehat{\mathcal{I}f}(0,0)-\hat{g}(0,0)\right|}^2+\sum _{v\in Q}{A}_v,$$

$$\begin{array}{rl}{A}_v& =\sum _{p\in \mathbb{Z}\setminus 0}{⟨pv⟩}^{2r}{\left|\widehat{\mathcal{I}f}(p{v}^{\perp },v)-\hat{g}(p{v}^{\perp },v)\right|}^2\\ & +\sum _{w\in {\mathbb{Z}}^2\setminus v\mathbb{Z}}{⟨w⟩}^{2r}{\left|\widehat{\mathcal{I}f}(w,v)-\hat{g}(w,v)\right|}^2.\end{array}$$

$$\begin{array}{rl}& {\left|\widehat{\mathcal{I}f}(0,0)-\hat{g}(0,0)\right|}^2+\alpha {\left|\hat{f}(0)\right|}^2\\ & +\sum _{v\in Q}\sum _{p\in \mathbb{Z}\setminus 0}\left({⟨pv⟩}^{2r}{\left|\hat{f}(p{v}^{\perp })-\hat{g}(p{v}^{\perp },v)\right|}^2+\alpha {⟨pv⟩}^{2s}{\left|\hat{f}(p{v}^{\perp })\right|}^2\right).\end{array}$$

$$\sum _{k\in {\mathbb{Z}}^2}{⟨k⟩}^{2r}\left({\left|\hat{f}(k)-\hat{g}(k,{\hat{k}}^{\perp })\right|}^2+\alpha {⟨k⟩}^{2(s-r)}{\left|\hat{f}(k)\right|}^2\right),$$

$$\hat{f}(k)=(1+\alpha {⟨k⟩}^{2(s-r)}{)}^{-1}\hat{g}(k,{\hat{k}}^{\perp }).$$

$$\underset{ϵ\to 0}{\lim }\underset{y\in Y}{\sup }\left\{\Vert {\mathcal{R}}_{\alpha (ϵ)}y-x{\Vert }_X;\Vert y-F(x){\Vert }_Y\le ϵ\right\}=0$$

Full text

This is the accepted manuscript of the work first published in SIAM J. Appl. Math., 80 (4), 1947–1976. 10.1137/19M1268070.

1. Introduction

We present a new computed tomography (CT) method for X-ray tomography in 2D, based on a relatively recent theoretical approach to X-ray tomography. The method reconstructs the Fourier series of a phantom via the projection of X-ray data into X-ray data on the flat torus which has a remarkably simple inverse X-ray transform. Therefore we call the new method Torus CT.

We emphasize that Torus CT is new as a numerical approach, so it is beyond the scope of a single article to give a full and detailed answer to all the natural questions that arise. This article should be seen as a proof of concept, exploring a new approach and showing that it is indeed feasible.

We have developed new mathematical theory and computational implementations.

We found a new inversion method and proved that it provides a regularization strategy with respect to suitable Sobolev norms; see sections 1.3. The numerical implementation was used to demonstrate the potential of Torus CT method in various simulations and tests, including data simulation in torus geometry and traditional experimental projections. Torus CT provided an efficient basis for inverse solution and its efficacy is shown in this work.

The article is organized as follows. In section 1.1 we give an overview of computed tomography and regularization, in section 1.2 we discuss works related to X-ray tomography on torus, and in section 1.3 we state the main theoretical results in this paper. Section 2 includes mathematical preliminaries, proofs of theorems and numerical analysis for Torus CT method. Section 3 contains mathematical formulation of computational forward and inverse models. Section 4 presents numerical experiments and their analysis. Conclusions are given in section 5. We have included a short note on supplementary material in the end of the article.

1.1. Overview of computed tomography and regularization methods

We give here an overview of X-ray tomography. Practical CT imaging was first introduced by Cormack and Hounsfield in 1970s based on the theoretical work of Cormark [3, 4] in 1960s. The mathematical theory itself was in fact earlier studied by Radon [28] in 1917. We give here only a narrow sample of topics and references in X-ray computed tomography. More references can be found in the cited works.

CT has many applications in medical imaging and engineering utilizing computerized axial tomography (CAT), positron-emission tomography (PET) and single-photon emission computed tomography (SPECT) [21]. Possible applications include imaging of patients in medicine and nondestructive testing in engineering. The most common inversion method for CT imaging is based on the filtered back-projection (FBP) algorithms [23, 14].

The FBP algorithms work well if there is sufficiently dense set of measurements. The filter in the standard FBP algorithm stands for applying the fractional Laplacian $(-\Delta)^{1/2}$ to the back-projected Radon transform data, and this does not include any regularization [9, 23]. In addition, regularization is required to control the errors in reconstructions caused by a measurement noise. See for example [22, 21].

Usually a regularization method is applied for a discretized X-ray tomography model as in the examples listed next. The most common regularization methods include Tikhonov regularization and truncated singular value decomposition (TSVD) which promote smoothness of reconstructions [22]. Other common regularization approaches include total variation (TV) regularization which promotes $\ell^{1}$ sparsity and, therefore, piecewise constant reconstructions [33, 8, 24, 7]. Another approach is to encode a priori information as a probability distribution and think the reconstruction problem as an Bayesian inverse problem for finding a posterior distribution [33, 17, 13, 8, 7].

The main difference of our proposed Tikhonov regularization approach, stated in theorem 2, compared to the usual regularization methods is that we do not discretize a phantom and regularization takes a form of a simple low-pass filter on the Fourier side. This also reflects the fact that Torus CT method is meshless (or meshfree) method. Theorem 3 states that the proposed regularization method is an admissible regularization strategy. Details are given in the subsequent sections.

1.2. The X-ray transform on torus, the Radon transform and the geodesic X-ray transform

In this paper we consider application of the X-ray transform on the flat torus $\mathbb{T}^{n}=\mathbb{R}^{n}/\mathbb{Z}^{n}$ to the usual CT in the case when $n=2$. In this section we give an account of theoretical works on the X-ray transforms on tori. As expected, the $d$-plane Radon transform of a function $f$ on $\mathbb{T}^{n}$ encodes the integrals of $f$ over all periodic $d$-planes. The X-ray transform corresponds to the case when $d=1$ and is in fact the geodesic X-ray transform on $\mathbb{T}^{n}$ over closed geodesics. It is described in section 2.3 how the usual CT reconstruction on $\mathbb{R}^{2}$ can be reduced to a reconstruction on $\mathbb{T}^{2}$.

Injectivity, reconstruction and certain stability estimates of the $d$-plane Radon transform on $\mathbb{T}^{n}$ were proved for distributions by Ilmavirta in [10]. The first injectivity result for the geodesic X-ray transform on $\mathbb{T}^{2}$ was obtained by Strichartz in [34], and generalized to $\mathbb{T}^{n}$ by Abouelaz and Rouvière in [2] if the Fourier transform is $\ell^{1}(\mathbb{Z}^{n})$. Abouelaz proved uniqueness under the same assumption for the $d$-plane Radon transform in [1]. A more general view and references on the Radon transform and the geodesic X-ray transform are given in [32, 9, 25, 12].

1.3. Inversion formulas and Tikhonov regularization

We state here our main theorems regarding the X-ray transform on $\mathbb{T}^{2}$. We write the X-ray transform on $\mathbb{T}^{2}$ as $\mathcal{I}$ and the formula

[TABLE]

holds true for sufficiently regular functions $f\colon\mathbb{T}^{2}\to\mathbb{C}$, including $L^{1}(\mathbb{T}^{2})$. In our proofs, we subsequently apply the Fourier slice theorem of the X-ray transform on $\mathbb{T}^{2}$, stated in formula (10). This fundamental, but simple, property is proved in [10]. The exact definition of the X-ray transform of the periodic distributions on $\mathbb{T}^{2}$ is given in section 2.

Our first theorem gives new inversion formulas for the X-ray transform. We give two proofs of theorem 1 in section 2.2. The first one does not rely on the inversion formula of [10] whereas the second simpler proof does.

Theorem 1**.**

Suppose that $f\in L^{1}(\mathbb{T}^{2})$. Let $k\in\mathbb{Z}^{2}$. If $k,v\neq 0$ and $v\bot k$, then

[TABLE]

If $k=0$, then

[TABLE]

The function $f$ can be reconstructed by Fourier series (9) and formulas (2) and (3).

Let $Q$ denote the set of all integer directions; a more detailed description will be given in section 2.4. We consider a Tikhonov minimization problem: given some data $g\in H^{r}(\mathbb{T}^{2}\times Q)$, find

[TABLE]

Let us define the post-processing operator $P^{s}_{\alpha}$ to be the Fourier multiplier $(1+\alpha\left\langle k\right\rangle^{2s})^{-1}$ and denote by $\mathcal{I}^{*}$ the adjoint of $\mathcal{I}$. We have the following theorems on regularization. The proofs are given in sections 2.4 and 2.5 respectively.

Theorem 2**.**

Let $r\in\mathbb{R}$, $s\geq r$, and $\alpha>0$. Suppose $g\in H^{r}(\mathbb{T}^{2}\times Q)$. The unique minimizer $f$ of the minimization problem (4) corresponding to Tikhonov regularization is $f=P^{s-r}_{\alpha}\mathcal{I}^{*}g\in H^{2s-r}(\mathbb{T}^{2})\subset H^{r}(\mathbb{T}^{2})$.

Theorem 3**.**

Suppose $r,t,s,\delta\in\mathbb{R}$ are such that $2s+t\geq r$, $\delta\geq 0$, and $s>0$. We assume that $f\in H^{r+\delta}(\mathbb{T}^{2})$ and $g\in H^{t}(\mathbb{T}^{2}\times Q)$.

Then our regularized reconstruction operator $P^{s}_{\alpha}\mathcal{I}^{*}$ gives a regularization strategy in the sense that

[TABLE]

where $\alpha(\varepsilon)=\sqrt{\varepsilon}$.

Moreover, if $\left\|g\right\|_{H^{t}(\mathbb{T}^{2}\times Q)}\leq\varepsilon$, $0<\delta<2s$ and $0<\alpha\leq 2s/\delta-1$, we have

[TABLE]

where $C(x)=x(x^{-1}-1)^{1-x}$.

Remark 4*.*

If we choose the regularization parameter as $\alpha=\varepsilon^{\gamma}$, the optimal asymptotic rate of convergence is obtained when $\gamma=(1+\delta/2s)^{-1}$. Then the terms $\alpha^{\delta/2s}$ and $\varepsilon/\alpha$ are of equal order.

We have also studied mapping properties, the adjoint and the normal operator of $\mathcal{I}$ in propositions 10 and 11; these results are stated in section 2. For example, it turns out that $\mathcal{I}|_{H^{s}(\mathbb{T}^{2})}$ is unitary to its range for any $s\in\mathbb{R}$ (see proposition 11).

Acknowledgements

J.I. was supported by the Academy of Finland (decision 295853). O.K. was supported by Jane and Aatos Erkko Foundation, and Jenny and Antti Wihuri Foundation. J.R. was supported by the Academy of Finland (CoE in Inverse Problems Research in 2017 & CoE in Inverse Modelling and Imaging in 2018–2019). The authors are grateful to Sampsa Pursiainen for providing access to a computer hardware used in simulations, Mikko Salo for helpful discussions, and Martin Bright for a helpful MathOverflow discussion related to proposition 16. O.K. wishes to thank Jari Hyttinen for providing a good working environment for completing this research. J.R. wishes to thank Matti Lassas and Samuli Siltanen for providing a good working environment at the University of Helsinki during his one-year visit in 2018. We thank the anonymous referees for their helpful comments.

2. Torus CT method

In this section we will lay out the theory of the Torus CT method. The reconstruction method is based on the Fourier series and properties of the geodesic X-ray transform on $\mathbb{T}^{2}$. There is a natural projection operator from the X-ray transform data of a compactly supported function on the plane to the X-ray transform data on $\mathbb{T}^{2}$. This so called torus-projection operator plays the role of the back-projection operator. For more details on the geodesic X-ray transform on tori see [10] and [11, Chapter 3].

2.1. The geodesic X-ray transform on $\mathbb{T}^{2}$

We define the flat torus as the quotient $\mathbb{T}^{2}:=\mathbb{R}^{2}/\mathbb{Z}^{2}$ and denote the quotient mapping $[\cdot]\colon\mathbb{R}^{2}\to\mathbb{T}^{2}$. A function $f\colon\mathbb{T}^{2}\to\mathbb{C}$ can be equivalently thought as a periodic function on $\mathbb{R}^{2}$ via the quotient mapping $[\cdot]$. We may thus consider a function $f\colon\mathbb{T}^{2}\to\mathbb{C}$ as a periodic function on the whole $\mathbb{R}^{2}$.

On closed Riemannian manifolds one defines the geodesic X-ray transform as a collection of line integrals of a function over periodic geodesics. All geodesics of $\mathbb{T}^{2}$ are given by the parametrizations $\gamma_{x,v}(t):=[x+tv],(x,v)\in[0,1]^{2}\times\mathbb{R}^{2}\setminus 0$. The geodesic $\gamma_{x,v}$ is periodic with period $1$ (with respect to the parameter $t$) if and only if $(x,v)\in[0,1]^{2}\times(\mathbb{Z}^{2}\setminus 0)$ (see e.g. [11, Exercise 23]). In general, a geodesic is periodic on $\mathbb{T}^{2}$ if and only if its direction vector is a multiple of a rational vector.

We denote the space of test functions by $\mathcal{T}:=C^{\infty}(\mathbb{T}^{2})$ and the set of all mappings $X\to Y$ by $Y^{X}$. We define the (geodesic) X-ray transform on $\mathbb{T}^{2}$ as an operator $\mathcal{I}\colon\mathcal{T}\to\mathbb{R}^{\mathbb{T}^{2}\times(\mathbb{Z}^{2}\setminus 0)}$ by

[TABLE]

A simple calculation shows that that $f\mapsto\mathcal{I}f(\cdot,v)$ is a formally self-adjoint operator on $\mathcal{T}$ for any fixed $v\in\mathbb{Z}^{2}\setminus 0$. We denote the dual space of $\mathcal{T}$ by $\mathcal{T}^{\prime}$, i.e. the space of distributions. By formal self-adjointness of $\mathcal{I}$, we may define the X-ray transform on distributions $f\in\mathcal{T}^{\prime}$ by

[TABLE]

where $(\cdot,\cdot)$ is the duality pairing.

If $f\in\mathcal{T}^{\prime}$, then we denote the Fourier coefficients of $f$ as $\hat{f}(k):=f(e^{-2\pi ik\cdot x}),k\in\mathbb{Z}^{2}$, and the Fourier series

[TABLE]

converges in the sense of distributions. We are now ready to recall the inversion formula in [10]:

Theorem 5** (Fourier slice theorem on $\mathbb{T}^{2}$, Eq. (9) in [10]).**

If $f\in\mathcal{T}^{\prime}$, then

[TABLE]

Theorem 5 gives a constructive formula (10) for the inverse of the X-ray transform on $\mathbb{T}^{2}$.

2.2. Inversion formula for integrable functions

In this section we simplify formula (10) for functions in $L^{1}(\mathbb{T}^{2})$. It turns out that the dimension of the integral defining $\widehat{\mathcal{I}f}(k,v)$ can be decreased by one using a change of coordinates, which enables a computationally faster implementation.

Recall that $f\in\mathcal{T}^{\prime}$ is in $L^{1}(\mathbb{T}^{2})$ if there exists a function $\tilde{f}\in L^{1}(\mathbb{T}^{2})$ such that

[TABLE]

It holds that $L^{1}(\mathbb{T}^{2})\subset\mathcal{T}^{\prime}$. If $\mathcal{I}f(\cdot,v)\in L^{1}(\mathbb{T}^{2})$ for some $f\in\mathcal{T}^{\prime}$, then we simply denote that $\mathcal{I}_{v}f=\mathcal{I}f(\cdot,v)$.

We define a family of coordinates which will be used repeatedly in this subsection. Suppose that $v\in\mathbb{Z}^{2}\setminus 0$ and $v_{1},v_{2}\neq 0$. Let $w_{m}:=\frac{m}{\left|v_{1}\right|}v,m\in\mathbb{Z}$, and define the coordinates $\varphi_{v,m}$ on $\mathbb{T}^{2}$ as

[TABLE]

Notice that $a\frac{v}{\left|v_{2}\right|}+(0,b)=(a\frac{v_{1}}{\left|v_{2}\right|},\frac{v_{2}}{\left|v_{2}\right|}(a+b))$ and $w_{m}=m(\frac{v_{1}}{\left|v_{1}\right|},\frac{v_{2}}{\left|v_{1}\right|})$ in the Cartesian coordinates. It easily follows that the Lebesgue measure on $\mathbb{T}^{2}$ transforms as $dm=\left|\frac{v_{1}}{v_{2}}\right|d(a,b)$ where $d(a,b)$ denotes the Lebesgue measure on $X:=[0,\left|\frac{v_{2}}{v_{1}}\right|)\times[0,1)$.

Remark 6*.*

The coordinates $\varphi_{v,m}$ parametrize $\mathbb{T}^{2}$ as parallelograms which are located on $\mathbb{R}^{2}$. Moreover, the parallelograms associated with $\varphi_{v,m},{m\in\mathbb{Z}}$ are disjoint for a fixed $v\in\mathbb{Z}^{2}\setminus 0$ when looked on $\mathbb{R}^{2}$. An example is given on Figure 1.

The next lemma states that the geodesic X-ray transform of an  $L^{1}(\mathbb{T}^{2})$ function can be defined geodesic-wise for almost every closed geodesic. Furthermore, the X-ray data for any fixed direction is also $L^{1}(\mathbb{T}^{2})$, and this definition agrees with the distributional definition.

Lemma 7**.**

Suppose that $v\in\mathbb{Z}^{2}\setminus 0$. Then the X-ray transform $\mathcal{I}_{v}\colon L^{1}(\mathbb{T}^{2})\to L^{1}(\mathbb{T}^{2})$ can be defined by the formula

[TABLE]

Moreover, we have:

1. (1)

This definition coincides with the distributional definition; for every $f\in L^{1}(\mathbb{T}^{2})$ and $g\in L^{\infty}(\mathbb{T}^{2})$ it holds that $(\mathcal{I}_{v}f,g)=(f,\mathcal{I}_{v}g)$. 2. (2)

$\mathcal{I}_{v}\colon L^{1}(\mathbb{T}^{2})\to L^{1}(\mathbb{T}^{2})$* is Lipschitz continuous with Lipschitz constant $1$.* 3. (3)

For almost every $p\in\mathbb{T}^{2}$ and every $v\in\mathbb{Z}^{2}\setminus 0$ and $t\in\mathbb{R}$ it holds that $\mathcal{I}_{v}f(p)=\mathcal{I}_{v}f(p+tv)$.

Proof.

This follows from Fubini’s theorem and straightforward calculations using the coordinates $\varphi_{v,m}$. We omit the details. ∎

We will give two proofs for theorem 1. The first proof is based on the assumption that $f\in L^{1}(\mathbb{T}^{2})$ and straightforward computation of the Fourier coefficients. The first proof proves the injectivity of the X-ray transform on $\mathbb{T}^{2}$ for $L^{1}(\mathbb{T}^{2})$ functions independently of [10]. The second proof is based on formula (10) and the assumption that $\mathcal{I}f(\cdot,v)\in L^{1}(\mathbb{T}^{2})$. Both of the proofs involve the coordinates $\varphi_{v,k}$.

First proof of theorem 1.

Recall that

[TABLE]

If $k_{1}=0$ or $k_{2}=0$, then formulas (2) and (3) follow trivially from (14).

The case $k_{1},k_{2}\neq 0$. We can use the coordinates $\varphi_{v,m},m\in\mathbb{Z}$, defined by formula (12). Using these coordinates we can calculate

[TABLE]

Notice that $k\cdot\left(\frac{a}{\left|v_{2}\right|}v+(0,b)+w_{m}\right)=k_{2}b$ since $v\cdot k=w_{m}\cdot k=0$.

Hence, we have

[TABLE]

We sum formula (16) for the values $m=0,\dots,\left|v_{1}\right|-1$, which gives

[TABLE]

This completes the proof. ∎

We will next prove a more general version of theorem 1.

Theorem 8**.**

Suppose that $f\in\mathcal{T}^{\prime}$ and $\mathcal{I}f(\cdot,v)\in L^{1}(\mathbb{T}^{2})$ for any $v\in\mathbb{Z}^{2}\setminus 0$. Then formulas (2) and (3) are true.

Proof.

We only show how to argue if $k_{1},k_{2}\neq 0$ since the other special cases are trivial. Recall that the inversion formula (10) states that $\widehat{\mathcal{I}_{v}f}(k)=\hat{f}(k)$ for any $v\in\mathbb{Z}^{2}\setminus 0$ such that $k\bot v$. We apply the coordinates $\varphi_{v,0}$.

Using Fubini’s theorem and calculations similar to the first proof of theorem 1, we get

[TABLE]

Now, formula (2) follows from property (3) of lemma 7. This completes the proof. ∎

Second proof of theorem 1.

Lemma 7 implies that $\mathcal{I}_{v}f\in L^{1}(\mathbb{T}^{2})$ if $f\in L^{1}(\mathbb{T}^{2})$. Hence, theorem 8 implies the inversion formulas. ∎

2.3. The torus-projection operator

We denote the X-ray transform of $f\colon\mathbb{R}^{2}\to\mathbb{C}$ by $\mathcal{R}_{v}f(p)$ for any $(p,v)\in\mathbb{R}^{2}\times S^{1}$. We parametrize the lines of the plane so that

[TABLE]

Suppose that $f$ is a compactly supported function on $\mathbb{R}^{2}$. We may then consider $f$ as a function defined on $\mathbb{T}^{2}$ after rescaling and periodizing. Let us denote the periodic extension of $f$ into $\mathbb{T}^{2}$ by the same symbol $f$.

Suppose further that $f\in C(\mathbb{T}^{2})$. As described in [11, Lemma 3.1], for any $p\in\mathbb{T}^{2}$ and $v\in\mathbb{Z}^{2}\setminus 0$ one can write $\mathcal{I}_{v}f(p)$ as a finite sum of terms $\mathcal{R}_{v/\left|v\right|}f(p_{i})$, $i=1,\dots,m$, and rescaling data by $\left|v\right|^{-1}$. One simply has to write any periodic geodesic $\gamma$ of $\mathbb{T}^{2}$ as a finite disjoint union of line segments that are supported in $[0,1)\times[0,1)$ and travel from the boundary to the boundary in the fundamental domain of $\mathbb{T}^{2}$. However, such unions are tedious to write down rigorously. This procedure defines the torus-projection operator $\mathcal{R}f\mapsto\mathcal{I}f$ for compactly supported continuous functions $f\colon\mathbb{R}^{2}\to\mathbb{C}$. We give example on Figure 2 which illustrates this procedure. For further details, see [11, Chapter 3]. See also the description of our numerical implementation in section 3.1.2. We have not studied mapping properties of the torus-projection operator, but we would like to see it studied thoroughly in the future. This would give new and important information about the connection between Euclidean and periodic X-ray transforms.

2.4. Sobolev spaces, adjoint, normal operator and regularization

Let $Q\subset\mathbb{Z}^{2}$ be such that every nonzero $v\in\mathbb{Z}^{2}$ is an integer multiple of a unique element in $Q$. We can simply take $Q$ to be the set of those vectors $(a,b)$ for which $a$ and $b$ are coprime with $a>0$ and $b\neq 0$ and the vectors $(0,1)$ and $(1,0)$. The set $Q$ is the set of all periodic directions on the torus, with all multiple counts removed. This set can be naturally identified with the projective space $\mathbb{P}^{1}$ defined later.

The X-ray transform we study takes a function on $\mathbb{T}^{2}$ to a function on $\mathbb{T}^{2}\times Q$. To set things up properly, we need to define function spaces and norms on both sides. On $\mathbb{T}^{2}$, we use the standard Sobolev scale of spaces $H^{s}(\mathbb{T}^{2})$ with the norms

[TABLE]

where $\left\langle k\right\rangle=(1+\left|k\right|^{2})^{1/2}$ as usual. On $\mathbb{T}^{2}\times Q$, we define the spaces $H^{s}(\mathbb{T}^{2}\times Q)$ to be the set of functions $g\colon\mathbb{T}^{2}\times Q\to\mathbb{C}$ for which

1. (i)

$g(\cdot,v)\in H^{s}(\mathbb{T}^{2})$ for every $v\in Q$; 2. (ii)

the average of every $g(\cdot,v)$ over $\mathbb{T}^{2}$ is the same; and 3. (iii)

the norm

[TABLE]

is finite. We set $v=0$ for the Fourier term $k=0$ to emphasize that it is the same for every $v\in Q$. We remind the reader that $0\notin Q$.

We emphasize that the regularity parameter $s$ can be any real number in the theory presented in this section. By setting $s=0$ one obtains a theory in $H^{0}=L^{2}$. We point out that the spaces considered here are different from [10].

Remark 9*.*

The norm of $g\in H^{s}(\mathbb{T}^{2}\times Q)$ is essentially an $\ell^{2}(Q)$ norm of the $H^{s}(\mathbb{T}^{2})$ norms of the functions $g(\cdot,v)$. This $\ell^{2}$ can be replaced with any $\ell^{p}$ without much effect to the theory, as the different functions in the family indexed by $Q$ have disjointly supported Fourier series apart from the origin. The case $p=\infty$ is particularly convenient because then special considerations are not needed at $k=0$. We choose $p=2$ to stay in a Hilbert space setting.

We denote $v^{\perp}=(-v_{2},v_{1})$ for any $v=(v_{1},v_{2})\in\mathbb{Z}^{2}$. For $v\in\mathbb{Z}^{2}\setminus 0$, we denote by $\hat{v}$ the unique point in $Q$ that is parallel to $v$. We can define $\hat{0}$ to be any point in $Q$; this choice will not matter. To keep notation neater, we will write $\hat{v}^{\perp}$ instead of $\widehat{v^{\perp}}$.

Proposition 10**.**

The X-ray transform is continuous $H^{s}(\mathbb{T}^{2})\to H^{s}(\mathbb{T}^{2}\times Q)$ for any $s\in\mathbb{R}$.

Proof.

For any $v\in Q$, the Fourier transform of function $\mathcal{I}f(\cdot,v)$ is supported on the line $v^{\perp}\mathbb{Z}$ by theorem 5. In fact, it is the restriction of $\hat{f}$ to this line. It then follows easily from the definition of the Sobolev norm on the Fourier side that $\mathcal{I}f(\cdot,v)\in H^{s}(\mathbb{T}^{2})$ whenever $f\in H^{s}(\mathbb{T}^{2})$.

It follows from the same theorem that $\widehat{\mathcal{I}f}(0,v)=\hat{f}(0)$ for all $v\in Q$, and so all the averages agree as required.

Since $\mathbb{Z}^{2}$ is a disjoint union of the origin and the punctured lines $v\mathbb{Z}\setminus 0$ with $v\in Q$, one can easily verify that $\left\|\mathcal{I}f\right\|_{H^{s}(\mathbb{T}^{2}\times Q)}=\left\|f\right\|_{H^{s}(\mathbb{T}^{2})}$. ∎

Proposition 11**.**

Fix any $s\in\mathbb{R}$. The adjoint of $\mathcal{I}\colon H^{s}(\mathbb{T}^{2})\to H^{s}(\mathbb{T}^{2}\times Q)$ is $\mathcal{I}^{*}\colon H^{s}(\mathbb{T}^{2}\times Q)\to H^{s}(\mathbb{T}^{2})$ given by

[TABLE]

The normal operator $\mathcal{I}^{*}\mathcal{I}\colon H^{s}(\mathbb{T}^{2})\to H^{s}(\mathbb{T}^{2})$ is the identity, so that $\mathcal{I}$ is unitary to its range.

Remark 12*.*

We emphasize that there is a striking difference with the usual Euclidean X-ray transform, where the normal operator is a convolution. In our setup the X-ray transform is directly inverted by its adjoint operator without any filtering or post-processing.

Proof of proposition 11.

Let us take any two functions $f\in H^{s}(\mathbb{T}^{2})$ and $g\in H^{s}(\mathbb{T}^{2}\times Q)$. We denote the complex conjugate of $z\in\mathbb{C}$ as $z^{*}$. Theorem 5 shows that $\widehat{\mathcal{I}f}(k,v)=\hat{f}(k)\delta_{0,k\cdot v}$, and so the $H^{s}$ inner products satisfy

[TABLE]

Therefore, the operator $\mathcal{I}^{*}$ defined above is the adjoint of $\mathcal{I}$.

It follows directly from the formula of theorem 5 that $\mathcal{I}^{*}$ is a left inverse of $\mathcal{I}$. ∎

Remark 13*.*

The X-ray transform or its normal operator have no effect on regularity. In the usual formulation, the normal operator does increase the smoothness index $s$, but when everything is set up on $\mathbb{T}^{2}$ the operators leave the regularity level intact.

We now turn to regularized inversion, and solve the Tikhonov minimization problem (4). We will make use of the post-processing operator $P^{s}_{\alpha}$, which is the Fourier multiplier $(1+\alpha\left\langle k\right\rangle^{2s})^{-1}$. It is evident that $P^{s}_{\alpha}$ maps continuously $H^{r}(\mathbb{T}^{2})\to H^{r+2s}(\mathbb{T}^{2})$ for any $s,r\in\mathbb{R}$.

Proof of theorem 2.

We begin with expanding the norms along lines given by $Q$ on the Fourier side. We have

[TABLE]

and

[TABLE]

where

[TABLE]

Each $\widehat{\mathcal{I}f}(w,v)$ vanishes in the last sum by theorem 5. Therefore, the second sum of $A_{v}$ is independent of $f$ and can be left out of the minimization problem. Furthermore, $\widehat{\mathcal{I}f}(pv^{\perp},v)=\hat{f}(pv^{\perp})$.

Thus, we are left with minimizing

[TABLE]

The notation introduced above allows us to rewrite the minimized quantity as

[TABLE]

and this can be minimized explicitly.

It suffices to choose each $\hat{f}(k)$ so that the term in the parentheses of (28) is minimized. A straightforward computation shows that the minimal $\hat{f}(k)$ is

[TABLE]

That is, the minimizer we sought is $f=P_{\alpha}^{s-r}\mathcal{I}^{*}g$. Finally, by the mapping properties of $P_{\alpha}^{s-r}$ and $\mathcal{I}^{*}$ we have $f\in H^{2s-r}(\mathbb{T}^{2})$. This implies that $f$ is in the correct space $H^{r}(\mathbb{T}^{2})$ since we assumed $s\geq r$. ∎

Remark 14*.*

Choosing $r=0$ and $s=1$, we reconstruct a function in $L^{2}$ with an $H^{1}$ penalty term. If we want the penalty to be the $L^{2}$ norm of the gradient without the $L^{2}$ norm of the function, the Fourier multiplier in the penalty term is changed from $\left\langle k\right\rangle^{2}$ to $\left|k\right|^{2}$. This corresponds to changing the Sobolev norm to a homogeneous Sobolev norm. Such changes lead to similar results but with slightly different postprocessing operator.

2.5. Regularization strategy

We define the concept of a regularization strategy according to [6, 15]. Let $X$ and $Y$ be subsets of Banach spaces and $F\colon X\to Y$ a continuous mapping. A family of continuous maps $\mathcal{R}_{\alpha}\colon Y\to X$ with $\alpha\in(0,\alpha_{0}]$ is called a regularization strategy if $\lim_{\alpha\to 0}\mathcal{R}_{\alpha}(F(x))=x$ for every $x\in X$. A choice of regularization parameter $\alpha(\epsilon)$ with $\lim_{\epsilon\to 0}\alpha(\epsilon)=0$ is called admissible if

[TABLE]

holds for every $x\in X$. Regularization strategies have been found for other inverse problems including, for example, electrical impedance tomography (EIT) [16] and inverse problem for the $1+1$ dimensional wave equation [18, 19].

We will next prove that the regularized inversion operator $P_{\alpha}^{s}\mathcal{I}^{*}$ obtained in theorem 2 actually provides an admissible regularization strategy with a quantitative stability estimate.

Proof of theorem 3.

Using proposition 11, we write

[TABLE]

and aim to estimate these two terms. In this proof, we denote the norm of $H^{r}(\mathbb{T}^{2})$ simply by $\left\|\cdot\right\|_{r}$.

Since $\left\|g\right\|_{H^{t}(\mathbb{T}^{2}\times Q)}\leq\varepsilon$ and $\left\|\mathcal{I}^{*}\right\|=\left\|\mathcal{I}\right\|=1$, we have $\left\|\mathcal{I}^{*}g\right\|_{t}\leq\varepsilon$. Applying the definitions of the norms and the operator $P^{s}_{\alpha}$, we find

[TABLE]

Estimating $1+\alpha\left\langle k\right\rangle^{2s}\geq\alpha\left\langle k\right\rangle^{2s}$ and using $-4s+2r-2t\leq 0$ shows that the supremum is at most $\alpha^{-2}$. Therefore

[TABLE]

which converges to zero as $\varepsilon\to 0$ with $\alpha=\sqrt{\varepsilon}$.

A calculation shows that $P_{\alpha}^{s}-\operatorname{id}=-\frac{\alpha\left\langle k\right\rangle^{2s}}{1+\alpha\left\langle k\right\rangle^{2s}}$ as a Fourier multiplier. Unfortunately, this implies that

[TABLE]

whenever $s>0$ and $\alpha>0$. Therefore, a uniform estimate is impossible when $\delta=0$, but it follows from the dominated convergence theorem that $\left\|(P^{s}_{\alpha}-\operatorname{id})f\right\|_{r}^{2}\to 0$ as $\alpha\to 0$ when $f\in H^{r}(\mathbb{T}^{2})$. The first claim of the theorem follows.

If $\delta>0$, the additional regularity of $f$ can be used to our advantage. It follows from the definitions of the norms that

[TABLE]

and thus

[TABLE]

Estimating this norm is crucial for the proof.

The supremum of (36) can be studied using the function $F\colon(0,\infty)\to(0,\infty)$ given by

[TABLE]

Simple calculus shows that if $2s>\delta$, then the maximum is attained at $x^{2s}=\alpha^{-1}(2s/\delta-1)$ and the maximal value on $(0,\infty)$ is

[TABLE]

We are interested in the maximum of $F$ on $[1,\infty)$. If $2s/\delta-1<\alpha$, then the maximum is reached at $x\in(0,1)$, and so the maximum on the relevant interval is $F(1)=\alpha/(1+\alpha)$. (One can also verify that the two maxima coincide when $2s/\delta-1=\alpha$, as they should.) We assumed that $2s/\delta>1$, so $\alpha\in(0,2s/\delta-1]$ for small enough $\alpha$.

For $\alpha\leq 2s/\delta-1$, the maximum value of $F$ is

[TABLE]

We conclude that

[TABLE]

and so

[TABLE]

Estimate (6) now follows easily from estimates (33) and (41). ∎

If $\alpha$ is bigger than assumed in the proof, then we may use the simpler estimate $F(x)\leq\alpha/(1+\alpha)\leq\alpha$ for all $x\geq 1$, which would lead to replacing $\alpha^{\delta/2s}C(\delta/2s)$ in estimate (6) by simply $\alpha$. However, we are only interested in the limit of small $\alpha$.

We point out that $C(\delta/2s)\to 1$ and $\alpha^{\delta/2s}\to 1$ when $\delta\to 0$, matching the norm in the limiting case of (34).

The noise $g$ in theorem 3 can be in any function space so that $\mathcal{I}^{*}g$, the reconstruction from pure noise, is in a suitable Sobolev space.

2.6. Numerical and asymptotic analysis for discretized problem

In this section, we consider questions arising from discrete practice. We analyze errors caused by a discretization of data in section 2.6.1. In section 2.6.2, we study how to choose a minimal set of X-ray directions in order to reconstruct all Fourier coefficients of a phantom in a given box.

Another source of errors in practice comes from the fact that we can only calculate finitely many coefficients of the Fourier series. The error caused by the cutoff of the Fourier series can be estimated with knowledge of asymptotic behavior of the Fourier coefficients. We do not consider this matter here further since it is a general question about convergence rates of Fourier series.

2.6.1. On convergence rates for discretization

Let $\mathbf{f}\in\mathbb{C}^{N}$ be written as $\mathbf{f}=(\mathbf{f}_{0},\dots,\mathbf{f}_{N-1})$. We define the discrete Fourier transform (DFT) of $\mathbf{f}$ by

[TABLE]

The following corollary of theorem 8 discretizes the inverse problem and reduces it to calculations of 1-dimensional DFTs. It is elementary and included here for completeness.

Corollary 15**.**

Let $f\in\mathcal{T}^{\prime}$, $\mathcal{I}_{v}f\in L^{1}(\mathbb{T}^{2})$, $k\in\mathbb{Z}^{2}\setminus 0$. Denote $g_{v}(y):=\mathcal{I}_{v}f(0,y)$ and $h_{v}(x):=\mathcal{I}_{v}f(x,0)$.

1. (1)

If $v\bot k$, then $\hat{f}(0,0)=\hat{g}_{(1,0)}(0)=\hat{h}_{(0,1)}(0)$ and

[TABLE] 2. (2)

(Left-point rule and DFT) Let $N\in\mathbb{Z}_{+}$. We denote $\textbf{g}_{l}=g_{v}(l/N)$ and $\textbf{h}_{l}=h_{v}(l/N)$ for $l=0,\dots,N-1$. If $\mathcal{I}_{v}f$ is Riemann integrable along vertical and horizontal lines, then

[TABLE]

Moreover, if $\mathcal{I}_{v}f\in C^{1}(\mathbb{T}^{2})$, then $\left|\hat{g}_{v}(k_{2})-\mathrm{DFT}(\textbf{g})_{k_{2}}\right|\leq C_{f,k_{2}}/N$ where $C_{f,k_{2}}>0$ does not depend on $N$. Similar statements hold for $h_{v}$ as well. 3. (3)

(Mid-point rule and DFT) Let $N\in\mathbb{Z}_{+}$. We denote $\textbf{g}_{l}=g_{v}(l/N+1/2N)$ and $\textbf{h}_{l}=h_{v}(l/N+1/2N)$ for $l=0,\dots,N-1$. If $\mathcal{I}_{v}f$ is Riemann integrable along vertical and horizontal lines, then

[TABLE]

Moreover, if $\mathcal{I}_{v}f\in C^{2}(\mathbb{T}^{2})$, then $\left|\hat{g}_{v}(k_{2})-\exp(-\pi ik_{2}/N)\mathrm{DFT}(\textbf{g})_{k_{2}}\right|\leq C_{f,k_{2}}/N^{2}$ where $C_{f,k_{2}}>0$ does not depend on $N$. Similar statements hold for $h_{v}$ as well.

Proof.

The statement (1) is a rephrased version of theorem 8. We only prove the statement (3). The proof of the statement (2) is similar and thus omited. Let $N\in\mathbb{Z}_{+}$ be fixed. By the definition of the DFT

[TABLE]

The statement follows since this the mid-point approximation of $\hat{g}_{v}(k_{2})$. The convergence rate is just a standard result on the mid-point rule (see e.g. [5]). ∎

2.6.2. Choosing directions for X-ray data.

Let us define the set

[TABLE]

where $\mathcal{Z}_{N}=[-N,N]^{2}\cap\mathbb{Z}^{2}$. Based on theorem 5, the data $(\mathcal{I}f(\cdot,v))_{v\in A_{N}}$ determines $(\hat{f}(k))_{k\in\mathcal{Z}_{N}}$. Thus, we define

[TABLE]

Define the set $V_{N}:=X_{+}\cup X_{-}\cup\{(1,0),(0,1)\}$ where

[TABLE]

Now, it is an elementary observation that the data $(\mathcal{I}f(\cdot,v))_{v\in V_{N}}$ determines $(\hat{f}(k))_{k\in\mathcal{Z}_{N}}$ and $\left|V_{N}\right|=\varphi(N)$.

We then turn to studying the asymptotic behavior of $\varphi(N)$. We denote by $\mathbb{P}^{1}:=\mathbb{P}^{1}(\mathbb{Q})$ the collection of equivalence classes $(a:b)$, $(a,b)\in\mathbb{Z}^{2}\setminus 0$, such that $(x,y)\in(a:b)$ if and only if $c(x,y)=(a,b)$ for some $c\neq 0$ and $(x,y)\in\mathbb{Z}^{2}\setminus 0$. The height is defined as $H(a:b):=\max\{\left|a\right|,\left|b\right|\}$ using the unique representative (up to a sign) of $(a:b)$ with $\gcd(a,b)=1$. One of the simplest special cases of Schanuel’s theorem [29, Theorem 1] states that

[TABLE]

as $N\to\infty$. More detailed exposition is given in the book of Serre [30, Chapter 2.5].

We conclude with the following proposition.

Proposition 16**.**

It holds that $\varphi(N)=\frac{2}{\zeta(2)}N^{2}+O(N\log N)$.

Proof.

If we want to reconstruct $\hat{f}(k)$, then we need at least one $v\in k^{\bot}$ by theorem 5 and, on the other hand, just one $v\in k^{\bot}$ is enough. It follows from the definition of height that

[TABLE]

The estimate follows now from Schanuel’s theorem. ∎

Remark 17*.*

The trivial estimate for directions needed in reconstruction of the Fourier coefficients $(\hat{f}(k))_{k\in\mathcal{Z}_{N}}$ would be $\varphi(N)\leq(2N+1)^{2}$. In comparison, proposition 16 implies that one needs to use asymptotically about $3/\pi^{2}\approx 30$ % of the data $(\mathcal{I}f(\cdot,v))_{v\in\mathcal{Z}_{N}}$.

3. Computational forward and inverse models

We have implemented two forward models for the X-ray transform on $\mathbb{T}^{2}$. The first forward model is based on direct integration over periodic geodesics on $\mathbb{T}^{2}$ (two different numerical integration schemes are implemented), and the second forward model on the usual Radon transform and the torus-projection operator. The regularized inverse model is based on theorems 1 and 2.

3.1. Computational forward models

3.1.1. Forward models on the torus

We have two different numerical integration schemes for the forward integration. The first one is analytical integration of a phantom which is discretized into square pixels of equal size. In this case, the forward operator, denoted by $\mathcal{A}_{1}$, is

[TABLE]

where $d_{i}$ is the length of the geodesic $\gamma_{x,v}$ and $f_{i}$ is the value of the discretized phantom in the $i$’th pixel, and $N$ is the size of the grid. The lengths $d_{i}$ are calculated by solving the intersection points of the line $\{\,x+tv\,;\,t\in[0,1]\,\}$ and the edges of the pixels when the pixels are periodically extended to $\mathbb{R}^{2}$.

In the second one, the integral is based on the use of global adaptive quadrature [31] which is implemented into Matlab’s integral function. In this case, a phantom is given in an analytical form. We denote this forward model by $\mathcal{A}_{2}$.

3.1.2. Forward model using the torus-projection and Radon data

This forward model corresponds to converting conventional X-ray data sets on $\mathbb{R}^{2}$ into X-ray data sets on $\mathbb{T}^{2}$. The forward model has two steps. The first step is to calculate Radon transform data using Matlab’s radon function. The second step is to calculate the torus-projection (see Section 2.3) of the Radon data. The directions for the Radon transform are chosen so that they contain all directions generated by integer vectors (see Section 2.6.2).

The X-ray beams on the radon function are parametrized by the distance between the line of a X-ray beam and the center of a domain $O$, and the angle of a X-ray beam measured from the $y$-axis into the counterclockwise direction. We denote simply that $Rf(v)=\text{{radon}}(f,\alpha_{v},M)$ where $\alpha_{v}$ is the angle defined above and $M$ is the number of X-rays taken into direction $v$. We index the rays as $k=1,\dots,M$. Further, denote the distances of rays to $O$ by $c_{k,v}$ and the projection values with the respective rays by $Rf(v)_{k}$.

We split each geodesic $\gamma_{x,v}$ into segments in which it travels from the boundary to the boundary when looked at the fundamental domain $[0,1]\times[0,1]$ of $\mathbb{T}^{2}$. Let $d_{x,v,i}$ be the distance of the $i$’th segment of the geodesic $\gamma_{x,v}$ and $O$, and $N$ the number of distinct segments. Finally, we can define the forward model $\mathcal{A}_{\mathbb{T}^{2}}$ as

[TABLE]

where

[TABLE]

if $\left|c_{k_{1,i},v}-d_{x,v,i}\right|+\left|c_{k_{2,i},v}-d_{x,v,i}\right|<\left|c_{k_{1,i},v}-c_{k_{2,i},v}\right|$, and $w_{1,i}=w_{2,i}=0$ otherwise.

The last condition ensures that the rays, corresponding to the data in interpolation, are on the different sides of the geodesic segment. Vice versa, if the condition does not hold, the geodesic segment is outside the projection width. In other words, this condition is the zero extension of the data near boundaries of the domain. In short, $\mathcal{A}_{\mathbb{T}^{2}}$ is the sum of weighted averages of two closest projection values with respect to their distances to the corresponding geodesic segments.

3.2. Computational inverse model

In the inverse model, we calculate the Fourier series coefficients of a phantom and reconstruct its Fourier series up to a finite radius $r>0$. The Fourier coefficients are calculated using the inversion formulas (2) and (3) of theorem 1. Furthermore, we Tikhonov regularize reconstructions using the filter $P_{\alpha}^{s}$ on the Fourier side according to theorems 2 and 3.

Let us write $B_{r}=B(0,r)\cap\mathbb{Z}^{2}$. The inverse model is

[TABLE]

where $\hat{f}_{\mathrm{rec}}(k)$ is calculated from data using the left-point rule and the DFT according to (1) and (2) of corollary 15. We remark that the inverse model is meshless and its output is a trigonometric polynomial.

4. Numerical experiments

4.1. Phantoms, convergence rates of Fourier series and discretization

4.1.1. Phantoms

We have used two phantoms in the numerical experiments, the Shepp–Logan phantom based on Matlab’s function phantom and the Flag phantom which is a piece-wise constant function representing a Nordic flag. The Flag phantom $f_{F}\colon[0,1]\times[0,1]\to\mathbb{R}$ was defined as

[TABLE]

where

[TABLE]

That is, $f_{F}$ describes the outer boundaries of the flag, and $g_{F}$ returns the background unless $x$ or $y$ is on the horizontal or vertical stripe, respectively.

4.1.2. Cutoff errors of Fourier series of phantoms

We analyzed the cutoff errors of Fourier series of the phantoms in order to determine a good, practical value of $r>0$ for the reconstructions. The squared cutoff error of Fourier series can be calculated via the formula

[TABLE]

using Parseval’s identity.

We computed $\epsilon_{r}$ for the Shepp–Logan phantom, the Flag phantom and the Flag phantom with a $45^{\circ}$ rotation. All the three phantoms were studied without noise and with salt-and-pepper (S&P) type noise applied to the phantoms using Matlab’s imnoise function with $0.02$ noise density. The phantoms were discretized into $4000\times 4000$ pixel grid and the Fourier coefficients $\hat{f}(k)$ were computed using Matlab’s fft2 and fftshift functions.

The squared cutoff errors $\epsilon_{r}$ are shown in Figure 3. The squared cutoff errors saturate at around $r=50$, though some improvement might be gained up to $r=200$. In our forward and inverse simulations, we have mainly used $r=50$ as it practically seems to be a sufficiently good choice.

4.1.3. Discretizations of phantoms and geodesics

The starting points of the used geodesics were chosen to be the equispaced points $\{(0,0),(1/n_{d},0),(2/n_{d},0),\dots,(1-1/n_{d},0)\}$ on the $x$-axis, except for geodesics in direction $v=(1,0)$ where the sampling was $\{(0,0),(0,1/n_{d}),(0,2/n_{d}),\dots,(0,1-1/n_{d})\}$ on the $y$-axis. In our experiments, we set $n_{d}=128$ when the cutoff radius of the Fourier series was $r\in\{50,100\}$, and $n_{d}=r$ when $r\in\{150,200\}$.

The phantoms were discretized with $512\times 512$ pixel grid when used for forward simulations with the forward models $\mathcal{A}_{1}$ and $\mathcal{A}_{\mathbb{T}^{2}}$. When we used the forward model $\mathcal{A}_{2}$, the Flag phantom was not discretized. The values of reconstructions were evaluated at equispaced points in $256\times 256$ pixel grid, and when compared to the ground truth, the Shepp–Logan and the Flag phantoms were discretized for the same grid.

4.2. Numerical analysis of forward models $\mathcal{A}_{1}$, $\mathcal{A}_{2}$ and $\mathcal{A}_{\mathbb{T}^{2}}$

4.2.1. Forward models $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ on the torus

We tested Torus CT using the Shepp–Logan phantom with simulated data

[TABLE]

We made reconstructions with cutoff radii $r\in\{50,100,150,200\}$ of the Fourier series.

In the case of $r=50$, we experimented with Tikhonov regularization. The reconstruction errors with different regularization parameters are shown in Figure 4. We have calculated the (relative) reconstruction errors using the formula

[TABLE]

The optimal regularization parameter values yielding the smallest error are given in Table 1. The plotted errors Figure 4 share some similarities in shape and the resulting regularization parameter values are close to each other.

The Shepp–Logan phantom is shown in Figure 5(a) and its non-regularized solution in Figure 5(b). The regularized solutions with $p=2$ and $p=\infty$ based regularization parameter values (Figures 5(d) and 5(e)) are similar, and $p=1$ based values yield slightly smoother reconstruction (Figure 5(c)).

We tested the effect of increasing the Fourier coefficient by computing the forward data required for reconstruction of the Fourier coefficients up to radii $r=100$, $r=150$ and $r=200$, and reconstructions are shown in Figures 5(f), 5(g), and 5(h) respectively. The constant regions in the phantom become a bit more smoother, but overall dynamical range is increased and the impact of noise in reconstructions remains relatively high.

Similar analysis was also performed with the Flag phantom. We simulated noisy data using the model $y=\mathcal{A}_{2}f+\mathcal{E}$ with the noise model of (55). The case $r=50$ was used to test regularization. The reconstruction errors $\epsilon_{p}^{\alpha,s}$ are shown in Figure 6 and the regularization values yielding the minimum error are given in Table 1. The regions close to the minimum of $\epsilon_{p}^{\alpha,s}$ are more distinct than in the case of the Shepp–Logan phantom, but similar shape is seen.

The Flag phantom is shown in Figure 7(a) and the non-regularized reconstruction in Figure 7(b). The regularized reconstructions with the optimal regularization parameters yielding the minimum errors with $p=1$, $p=2$ and $p=\infty$ are shown in Figures 7(c), 7(d) and 7(e), respectively. The regularization parameter values yielding the minimum were close to each other, and with the Flag phantom, no significant difference is seen in the regularized reconstructions.

Increasing the radius of the Fourier coefficients again increases the dynamical range, plotted in Figures 7(f), 7(g) and 7(h) for $r\in\{100,150,200\}$, respectively. However, unlike with the Shepp–Logan phantom, the details become more distinct, especially the details of the corners in the Flag phantom.

4.2.2. Forward model $\mathcal{A}_{\mathbb{T}^{2}}$ using the torus-projection and Radon data

To test how a Torus CT algorithm would work with experimental data acquisition, we computed Radon transform of the phantoms and projected it to $\mathbb{T}^{2}$ using the model $\mathcal{A}_{\mathbb{T}^{2}}$ with noise on each data point on $Rf(v)_{k}$. More precisely, we simulated data according to formula (50) where each $Rf(v)_{k}$ was replaced by noisy data $Rf(v)_{k}+\mathcal{E}$ where $\mathcal{E}\sim\mathcal{N}(0,\sigma^{2})$ with $\sigma=0.02$.

The projection directions for Radon transform were computed such that they determined the Fourier coefficients up to radius $r=50$. An illustration of how the projection directions in $(0,90)^{\circ}$ are distributed is shown in Figure 8, and the remaining projection directions are reflections of the projection directions in $(0,90)^{\circ}$ about the $y$-axis. In total, with $r=50$, there are 3097 unique projection directions. Two major concentrations of the directions are close to $45^{\circ}$, both above and below, but also smaller concentrations are found elsewhere, e.g., close to $22.5^{\circ}$.

The reconstructions from data computed with $\mathcal{A}_{\mathbb{T}^{2}}$ are shown in Figure 9. Shepp–Logan (Figures 9(a) and 9(d)) and $30^{\circ}$ rotated Flag (Figures 9(c) and 9(f)) are reconstructed well even with the noisy data. However, the non-rotated Flag phantom (Figures 9(b) and 9(e)) is rather poor.

As Figure 8 already showed, the torus-optimized data does not contain information uniformly of all directions. Instability and thus blurriness is expected where the available directions shown in Figure 8 are sparse. The most prominent example is the coordinate axis. Hence, it was expected that recovery of singularities near coordinate axis directions is practically unstable [26, 27]. This unfavorable feature has motivated us to study a modification of the inversion method which is expected to succeed better in practice. This modified approach and the related simulations are presented in section 4.3. With the Shepp–Logan phantom, the features are clearly detected, especially in the noise-free case (Figure 9(a)) indicating potential in the technique. Regularized solutions are shown in Figures 9(g)), 9(h) and 9(i) from the Shepp–Logan, the non-rotated and the rotated Flag phantoms, repectively. The regularization smoothed the reconstructions, decreased their dynamic range and no additional features were revealed from the noise. The regularization parameter values were $\alpha=0.75$ and $s=0.5$, chosen with manual experimentation.

For comparison, we computed the respective FBP reconstructions (shown in Figure 10) with Matlab’s iradon function using default settings. The projection data $R_{\texttt{radon}}f(v)_{k}+\mathcal{E}$ was down sampled by factor of $2$ with imresize to match reconstruction resolution $256\times 256$. It seems, that the uneven distribution of projection angles creates errors in reconstruction, since similar artefacts in horizontal, vertical and diagonal directions are seen also in the FBP reconstruction as in the ones computed with the Torus CT method in Figure 9. From the FBP this was expected as it is prone to streaking. In general, the FBP reconstructions are of good quality, since there is a lot of data available. With the same number of projections, 3097, but evenly distributed as they normally are, the FBP reconstruction are better quality than any other presented in this paper.

The error $\epsilon_{p}=\|f-f_{\mathrm{rec}}^{\mathrm{FBP}}\|_{L^{p}(\mathbb{R}^{2})}/\|f\|_{L^{p}(\mathbb{R}^{2})}$ between the FBP reconstruction $f_{\mathrm{rec}}^{\mathrm{FBP}}$ and the phantom $f$ is tabulated in Table 2. When compared with $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ and Shepp–Logan and Flag phantoms, with all values of $p$, the errors $\epsilon_{p}$ are higher than errors of regularized Torus CT reconstructions $\epsilon_{p}^{\alpha,s}$ shown in Table 1. When compared with the errors of non-regularized reconstructions $\epsilon_{1}^{0,0}$ and $\epsilon_{2}^{0,0}$, the FBP and Torus CT are relatively close, but the error $\epsilon_{\infty}^{0,0}$ of Torus CT is lower with the Shepp–Logan phantom and higher with the Flag phantom.

For use with practical data acquisition with commonly used equispaced projection angles, Torus CT requires more work to handle the increased additive noise in $\mathcal{A}_{\mathbb{T}^{2}}$, since the reconstructions have more noise outside of the support of the phantom than in the FBP reconstructions as seen in Figure 9. The respective errors $\epsilon_{1}^{\alpha,s}$ and $\epsilon_{2}^{\alpha,s}$ are higher than those of the FBP, presented in Table 2. Nevertheless, in terms of reconstructing the correct dynamical range of the objects, measured with $\epsilon_{\infty}^{\alpha,s}$ and $\epsilon_{\infty}$, the Torus CT method is equivalent or better than with the FBP.

4.3. Rotating phantom on torus

The theoretical formulation allows to place a phantom inside $\mathbb{T}^{2}$ in many different positions. Such choice of an orientation of a phantom leads to different choice of projection directions, and thus results different reconstructions when only finitely many Fourier coefficients are recovered. This motivated to test, whether the reconstructions improve when data is acquired from several rotated phantoms. We verified this hypothesis by computing the data of the Shepp–Logan phantom from nine different rotational orientations with $20^{\circ}$ interval and the Flag phantom with six different rotational orientations with $30^{\circ}$ interval, both with Fourier coefficient radius $r=50$.

Denote the angles of rotational data sets of the Shepp–Logan phantom with

[TABLE]

and of the Flag phantom with

[TABLE]

where $i\in\mathbb{N}$. Furthermore, denote rotations about the point $(0.5,0.5)$ by the symbols $\theta_{k}\in\Theta_{i}$.

The forward solutions were computed for the rotated phantoms as described in subsection 4.2. Then, we reconstructed the rotated phantoms using formula (51) and rotated the reconstructions back to the natural orientation, i.e., $-\theta_{i}$ about the point $(0.5,0.5)$. We denote these naturally oriented reconstructions by $f_{\text{rec}}^{\alpha,s}(x;k)$ for each $k=1,\dots,i$. Eventually, we computed the average of the rotated reconstructions

[TABLE]

For the Shepp–Logan phantom, the rotation was computed using Matlab’s imrotate with crop option for the Shepp–Logan phantom. For the Flag phantom, rotation about the point $(0.5,0.5)$ was made by changing the coordinates in equations (52) and (53). For both phantoms, the rotation of the reconstruction to natural orientation was computed using imrotate.

The reconstructions $f_{\Theta_{i}}$ are shown in Figure 11 for the Shepp–Logan phantom and in Figure 12 for the Flag phantom. The reconstruction errors $\epsilon_{p}^{0,0}$ (eq. 56) are tabulated in Table 3 for both phantoms. With the Flag phantom, the reconstruction errors are computed only on the support of the flag phantom and the errors are ruled out at the points where the Flag phantom vanish. With the Shepp–Logan, $\epsilon_{p}^{0,0}$ was computed on the whole grid. The reconstructions were evaluated in a grid of $256\times 256$ pixels.

With the Shepp–Logan phantom there is clear visual improvement with the increase of rotationally acquired data and also decreasing trend in the errors. With the simpler Flag phantom, reconstruction quality seems to saturate as the decrease of error norms is not as clear as with the Shepp–Logan phantom. Nevertheless, the smallest error norm is given with the highest number of rotational data (shown in Table 3) and the visual evaluations support this.

The rotation of the phantom does contribute to the improvement of reconstructions. In other words, the improvement is not merely due to averaging out the zero mean noise. We simulated data from the same rotational orientation of the Shepp–Logan phantom nine times, and the errors were as follows: $\epsilon_{1}^{0,0}=52\%$, $\epsilon_{2}^{0,0}=37\%$ and $\epsilon_{\infty}^{0,0}=66\%$. For the Flag phantom with six times from same rotational orientation, the errors were: $\epsilon_{1}^{0,0}=43\%$, $\epsilon_{2}^{0,0}=28\%$ and $\epsilon_{\infty}^{0,0}=76\%$. In both cases the errors were higher than what was gained with different rotational orientations, except for $\epsilon_{\infty}^{0,0}$ with Shepp–Logan phantom where the error was almost equal. It should also be noted that the use imrotate induces some blurring during the rotation of reconstructions and the Shepp–Logan phantom. Nonetheless, rotational reconstructions approach performed better.

4.4. Computing times

The computing time of the forward system, i.e., data, depends mainly on the cutoff radius $r$ of the Fourier series and on the number used geodesics in each direction (see discretization of $x$ in 4.1.3). In terms of this paper, the radius $r$ was more of the interest. Discretization relates to the numerical accuracy and the data acquisition accuracy of experimental setup. Example computing times $t_{r}$ for data on Lenovo P51 laptop with Intel i7-7820HQ CPU and 32 GB of RAM having MATLAB R2017a (The MathWorks, Inc.) with the Shepp–Logan phantom and $\mathcal{A}_{1}$ are $t_{50}=5.5\leavevmode\nobreak\ \text{min}$ and $t_{100}=51\leavevmode\nobreak\ \text{min}$; and with the Flag phantom and $\mathcal{A}_{2}$, $t_{50}=100\leavevmode\nobreak\ \text{min}$ and $t_{100}=62\leavevmode\nobreak\ \text{min}$. On Lenovo P910 high-end workstation with two Intel Xeon E5-2697 processors and 256 GB RAM having MATLAB version R2016b 64-bit (The MathWorks, Inc.), the computation times were of order $t_{50}=2.2\leavevmode\nobreak\ \text{min}$, $t_{100}=15\leavevmode\nobreak\ \text{min}$, $t_{150}=60\leavevmode\nobreak\ \text{min}$ and $t_{200}=188\leavevmode\nobreak\ \text{min}$ with the Shepp–Logan phantom and $\mathcal{A}_{1}$; and $t_{50}=1.6\leavevmode\nobreak\ \text{min}$, $t_{100}=21\leavevmode\nobreak\ \text{min}$, $t_{150}=79\leavevmode\nobreak\ \text{min}$ and $t_{200}=242\leavevmode\nobreak\ \text{min}$ with the Flag phantom $\mathcal{A}_{2}$. The analytical integration applied when using the Shepp–Logan phantom with $\mathcal{A}_{1}$ explains its faster computations times.

The projection of Radon transform sinogram to the torus $\mathcal{A}_{\mathbb{T}^{2}}$ and its reconstruction lasted approximately eight minutes on Lenovo P51. However, the current implementation was not optimized at all and included, among other, three nested for-loops. Hence, here the computational efficiency will increase during further development.

5. Conclusions

Our main conclusion is that the new inversion method for the X-ray transform based on the torus is numerically applicable. We emphasize, however, that this is only a first step. Much is left to be studied, such as the details of the mapping that transforms the planar X-ray transforms onto the torus and application of different regularization methods. We have developed the previously unapplied theory in a direction relevant to computation and carried out first numerical tests.

Our theoretical results were strongly motivated by practical requirements, including new and computationally fast reconstruction formulas from X-ray data in theorem 1, and rigorous mathematical theory for Tikhonov regularized reconstructions from X-ray data on the flat torus in theorems 2 and 3. We gave mathematical formulations of discretized forward and inverse models in section 3, and considered numerical analysis in section 2.6. The numerical implementation using Matlab was described in section 4.

The numerical implementation demonstrated the efficacy of Torus CT. Torus CT is computationally relatively efficient compared with iterative techniques, though still slower than current implementations of the FBP. An interesting feature of Torus CT is its meshless nature: Once the Fourier coefficients are computed, the reconstruction can be evaluated in any desired grid points. Currently, theory and the implementation are established in 2 dimensions, which is suitable for slice-wise reconstructions of 3-dimensional structures. One future research direction could be the development of algorithms and theory in higher dimensional settings.

Data simulation was also computed with the traditional Radon transform corresponding to experimental image acquisition with projection angles preferred by the Torus CT. Reconstruction quality was promising but not quite competitive with state-of-the-art methods. Some initial work has been conducted with conventional, evenly distributed projection angles, in which case there are various ways to interpolate projection data to directions preferred by Torus CT. It seems that rotations of a phantom could result sharp reconstructions and might allow reduction in the number of projection directions. This question should be studied more and with experimental X-ray data.

We find that Torus CT opens up a promising path to studying X-ray tomography in a new way. It provides a useful setup for future study — both theoretical, numerical, and practical — of the many remaining questions.

Supplementary material

Matlab code

We provide Creative Commons 4.0 licensed Matlab code files that implement forward model $\mathcal{A}_{1}$ (section 3.1.1) and inverse solutions on torus (section 3.2). The code package comprises of three files: TorusCTrun.m is the main script, DFT.m implements discrete Fourier transform (section 2.6), and LineIntegralOnGrid.m computes the exact line integral (49) over periodically extended, pixelized phantom. The files are available at [20].