# Ellipsephic harmonic series revisited

Jean-Paul Allouche
CNRS, IMJ-PRG, Sorbonne, 4 Place Jussieu
F-75252 Paris Cedex 05, France
<EMAIL_ADDRESS>
Yining Hu
Institute for Advanced Study in Mathematics
Harbin Institute of Technology
Harbin, PR China
<EMAIL_ADDRESS>
Claude Morin
<EMAIL_ADDRESS>
###### Abstract

Ellipsephic or Kempner-like harmonic series are series of inverses of integers
whose expansion in base $B$, for some $B\geq 2$, contains no occurrence of
some fixed digit or some fixed block of digits. A prototypical example was
proposed by Kempner in 1914, namely the sum inverses of integers whose
expansion in base $10$ contains no occurrence of a nonzero given digit.
Results about such series address their convergence as well as closed
expressions for their sums (or approximations thereof). Another direction of
research is the study of sums of inverses of integers that contain only a
given finite number, say $k$, of some digit or some block of digits, and the
limits of such sums when $k$ goes to infinity. Generalizing partial results in
the literature, we give a complete result for any digit or block of digits in
any base.

Keywords: Kempner-like series. Ellipsephic numbers. Sum of digits. Counting
blocks of digits.

MSC Classification: 11A63, 11B85, 68R15.

##  1 Introduction

While the harmonic series $\sum\frac{1}{n}$ is divergent, restricting the
indices in the sum to integers satisfying innocent-looking conditions can
yield convergent series. One of the first such examples is probably the 1914
result of Kempner [21] stating that the sum of inverses of integers whose
expansion in base $10$ contains no occurrence of a given digit ($\neq 0$)
converges. After Kempner’s paper and the 1916 paper of Irwin [20], several
papers addressed extensions or generalizations of this result, as well as
closed forms of numerical computations of the sums of the corresponding
series: see, e.g., [2, 9, 10, 11, 12, 14, 15, 16, 17, 22, 23, 24, 25, 26, 27,
28, 29, 30, 31, 32, 33] and the references therein.

We will revisit harmonic series with missing digits or missing blocks of
digits: these series are called Kempner-like harmonic series or ellipsephic
harmonic series in the literature (for the origin of the term ellipsephic
coined by C. Mauduit, one can look at [13, p. 12] and [18, Footnote, p. 6];
also see the discussion in [6]). More precisely, let $B$ an integer $\geq 2$
and let $a_{w,B}(n)$ denote the number of occurrences of the word (string) $w$
in the base $B$ expansion of the integer $n$ (where $|w|$ is the length of the
string $w$), and $s_{B}(n)$ the sum of the digits of the base $B$ expansion of
$n$. It was proven by Farhi [15] that, for any digit $j$ in
$\\{0,1,\dots,9\\}$, one has

$\lim_{k\to\infty}\sum_{\begin{subarray}{c}n\geq 1\\\
a_{j,10}(n)=k\end{subarray}}\frac{1}{n}=10\log 10.$

As explained in [6], a post on mathfun asked for the value of

$\displaystyle\lim_{k\to\infty}\sum_{\begin{subarray}{c}n\geq 1\\\
s_{2}(n)=k\end{subarray}}\frac{1}{n}\cdot$

It was proved in [6] that

$\lim_{k\to\infty}\sum_{\begin{subarray}{c}n\geq 1\\\
s_{B}(n)=k\end{subarray}}\frac{1}{n}=\frac{2\log B}{B-1}\cdot$ (1.1)

Thus we have of course

$\lim_{k\to\infty}\sum_{\begin{subarray}{c}n\geq 1\\\
a_{1,2}(n)=k\end{subarray}}\frac{1}{n}=\lim_{k\to\infty}\sum_{\begin{subarray}{c}n\geq
1\\\ s_{2}(n)=k\end{subarray}}\frac{1}{n}=2\log 2.$

Here we will evaluate all such series, where $a_{1,B}(n)$ is replaced with
$a_{w,B}(n)$, where $w$ is any block of digits in base $B$, by proving that

$\lim_{k\to\infty}\sum_{\begin{subarray}{c}n\geq 1\\\
a_{w,B}(n)=k\end{subarray}}\frac{1}{n}=B^{|w|}\log B.$ (1.2)

To prove this result, we will replace $1/n$ with the seemingly more
complicated function $\log_{B}(\frac{n}{n+1})$, and make use of results proved
in (or inspired by) [4, 19]. In passing we will generalize [4, 6] and re-prove
[19].

##  2 “Reducing” the problem

Let $B$ be an integer $\geq 2$. Let $w$ be a string of letters in
$\\{0,1,\dots,B-1\\}$. Let $a_{w,B}(n)$ be the number of (possibly
overlapping) occurrences of of the string $w$ in the base $B$ expansion of
$n$. First we note that the series
$\displaystyle\sum_{\begin{subarray}{c}n\geq 1\\\
a_{w,B}(n)=k\end{subarray}}\frac{1}{n}$ converges: the proof is the same as in
[4, Lemma 1, p. 194], namely one uses a counting argument for the case of a
single digit, and one replaces the base with some of its powers for the case
of a block of digits. Now, to evaluate the series, the idea is to replace it
with a convergent series $\displaystyle\sum_{\begin{subarray}{c}n\geq 1\\\
a_{w,B}(n)=k\end{subarray}}b_{w}(n)$ whose sum, say $A_{w}(k)$, tends to a
limit, say $A_{w}$ when $k\to\infty$. Furthermore, if we have the property
$b_{w}(n)-1/(B^{|w|}n)={\mathcal{O}}_{w}(1/n^{2})$ when $n$ tends to infinity,
then we obtain

  * $\sum_{\begin{subarray}{c}n\geq 1\\\ a_{w,B}(n)=k\end{subarray}}\frac{1}{n}\ \text{converges, and \ \ }$

  * $\lim_{k\to\infty}\sum_{\begin{subarray}{c}n\geq 1\\\ a_{w,B}(n)=k\end{subarray}}\frac{1}{n}=B^{|w|}A_{w}.$

(Note that if $a_{w,B}(n)$ tends to infinity, then $n$ must also tend to
infinity.)

Inspired by [4, 19], we define $L(n)$ by $L(0):=0$ and
$L(n):=\log_{B}\left(\dfrac{n}{n+1}\right)$ for $n\geq 1$. For a string $w$
over the alphabet $[0,B-1]$, let $v(w)$ denote the integer whose expansion in
base $B$ is $w$ (with possible leading $0$’s if $w\neq 0^{j}$).

######  Proposition 1.

Let $w$ be a nonempty string over the alphabet $[0,B-1]$,

$g=B^{|w|-1},\quad h=\left\lfloor\frac{v(w)}{B}\right\rfloor.$

Then, for all $k\geq 0$,

$\sum_{\begin{subarray}{c}n\\\ a_{w,B}(gn+h)=k\end{subarray}}L(Bgn+v(w))=-1,$
(2.1)

where the sum is over $n\geq 1$ if $w=0^{j}$ and $n\geq 0$ otherwise.

###### Proof.

Let $c$ be the last letter of $w$. Let $d_{w}(k)$ be defined by

$d_{w}(k)=\sum_{\begin{subarray}{c}n\geq 0\\\
a_{w,B}(n)=k\end{subarray}}L(Bn+c).$

(Note that this series converges since $L(n)\sim\frac{1}{n\log B}$ when $n$
goes to infinity.) By writing $n=gr+m$, with $r\geq 0$ and $0\leq m\leq g-1$,
we see that

$d_{w}(k)=\sum_{m=0}^{g-1}\sum_{\begin{subarray}{c}r\geq 0\\\
a_{w,B}(gr+m)=k\end{subarray}}L(Bgr+Bm+c).$

Similarly, if we let

$e_{w}(k)=\sum_{\begin{subarray}{c}n\geq 0\\\
a_{w,B}(Bn+c)=k\end{subarray}}L(Bn+c)$

(which is convergent, like $d_{w}(k)$), then we have

$e_{w}(k)=\sum_{m=0}^{g-1}\sum_{\begin{subarray}{c}r\geq 0\\\
a_{w,B}(Bgr+Bm+c)=k\end{subarray}}L(Bgr+Bm+c).$

Note that

$a_{w,B}(Bgr+Bm+c)-a_{w,B}(gr+m)=\begin{cases}1&\mbox{if }m=h\\\
0&\mbox{otherwise}\end{cases}$

for $r\geq 0$ if $w\neq 0^{j}$ and for $r\geq 1$ if $w=0^{j}$. For $w=0^{j}$
and $r=0$, the above difference is $0$ for all $m$ because we do not pad
leading $0$’s in this case.

Therefore

$d_{w}(k)-e_{w}(k)=\sum_{\begin{subarray}{c}r\\\
a_{w,B}(gr+h)=k\end{subarray}}L(Bgr+v(w))\\\ -\sum_{\begin{subarray}{c}r\\\
a_{w,B}(gr+h)=k-1\end{subarray}}L(Bgr+v(w))$ (2.2)

where the sum is over $r\geq 1$ if $w=0^{j}$ and $r\geq 0$ otherwise.

If we could show that $d_{w}(k)=e_{w}(k)$ for $k>0$, then it would follow from
equation (2.2) that the value of the sum

$\sum_{\begin{subarray}{c}r\\\ a_{w,B}(gr+h)=k\end{subarray}}L(Bgr+v(w))$

is independent of $k$ and hence equal to $d_{w}(0)-e_{w}(0)$. To prove this,
notice that

$L(n)-\sum_{j=0}^{B-1}L(Bn+j)=\begin{cases}0&\mbox{ if }n>0\\\ 1&\mbox{ if
}n=0,\end{cases}$ (2.3)

and

$\displaystyle\quad\sum_{\begin{subarray}{c}n\geq 0\\\
a_{w,B}(n)=k\end{subarray}}L(n)$
$\displaystyle=\sum_{j=0}^{B-1}\sum_{\begin{subarray}{c}n\geq 0\\\
a_{w,B}(Bn+j)=k\end{subarray}}L(Bn+j)$
$\displaystyle=\\!\\!\\!\sum_{\begin{subarray}{c}n\geq 0\\\
a_{w,B}(Bn+c)=k\end{subarray}}L(Bn+c)+\sum_{\begin{subarray}{c}j=0\\\ j\neq
c\end{subarray}}^{B-1}\sum_{\begin{subarray}{c}n\geq 0\\\
a_{w,B}(n)=k\end{subarray}}L(Bn+j).$

Hence

$\displaystyle e_{w}(k)$ $\displaystyle=\sum_{\begin{subarray}{c}n\geq 0\\\
a_{w,B}(Bn+c)=k\end{subarray}}L(Bn+c)$
$\displaystyle=\sum_{\begin{subarray}{c}n\geq 0\\\
a_{w,B}(n)=k\end{subarray}}(L(n)-\sum_{\begin{subarray}{c}j=0\\\ j\neq
c\end{subarray}}^{B-1}L(Bn+j))$ $\displaystyle=\sum_{\begin{subarray}{c}n\geq
0\\\
a_{w,B}(n)=k\end{subarray}}(L(n)-\sum_{\begin{subarray}{c}j=0\end{subarray}}^{B-1}L(Bn+j))+\sum_{\begin{subarray}{c}n\geq
0\\\ a_{w,B}(n)=k\end{subarray}}L(Bn+c)$

By (2.3), the first sum is $1$ if $k=0$, and $0$ if $k>0$. The second sum is
the definition of $d_{w}(k)$. ∎

######  Lemma 2.

Let $t$ be an integer whose expansion in base $B$ is $t=b_{1}b_{2}\ldots
b_{s}$. For $1\leq r\leq s$,

If $b_{1}\ldots b_{r}$ is not a suffix of $w$, then

$\sum_{\begin{subarray}{c}n\\\ a_{w,B}(B^{r}n+v(b_{1}\ldots
b_{r}))=k\end{subarray}}L(B^{s}n+t)=\sum_{\begin{subarray}{c}n\\\
a_{w,B}(B^{r-1}n+v(b_{1}\ldots b_{r-1}))=k\end{subarray}}L(B^{s}n+t).$

If $b_{1}\ldots b_{r}$ is a suffix of $w$, then

$\sum_{\begin{subarray}{c}n\\\ a_{w,B}(B^{r}n+v(b_{1}\ldots
b_{r}))=k\end{subarray}}L(B^{s}n+t)=\\\ \sum_{\begin{subarray}{c}n\\\
a_{w,B}(B^{r-1}n+v(b_{2}\ldots
b_{r}))=k\end{subarray}}L(B^{s-1}n+t^{\prime})\\\
-\sum_{\begin{subarray}{c}j=0\\\ j\neq
b_{1}\end{subarray}}^{B-1}\sum_{\begin{subarray}{c}n\\\
a_{w,B}(B^{r-1}n+v_{j})=k\end{subarray}}L(B^{s}n+t_{j})$ (2.4)

where $t^{\prime}=v(b_{2}\cdots b_{s})$, $t_{j}=v(jb_{2}\cdots b_{s})$, and
$v_{j}=v(jb_{2}\ldots b_{r-1})$ if $r\geq 2$ and $v_{j}=0$ if $r=1$.

The proof is the same as in [4].

######  Theorem 3.

There is a rational function $b_{w}(n)$ such that for all $k\geq 0$ we have

$\sum_{\begin{subarray}{c}n\\\
a_{w,B}(n)=k\end{subarray}}\log_{B}(b_{w}(n))=-1.$ (2.5)

(The summation is over $n\geq 1$ for $w=0^{j}$ and $n\geq 0$ otherwise.) and

$\log(b_{w}(n))=-\frac{1}{B^{|w|}n}+\mathcal{O}(1/n^{2}).$ (2.6)

###### Proof.

The existence of $b_{w}(n)$ follows from Proposition 1 and iterated
applications of Lemma 2. The process of obtaining $b_{w}(n)$ can be visualized
by a tree $T$ whose root is

$\sum_{\begin{subarray}{c}n\\\ a_{w,B}(gn+h)=k\end{subarray}}L(Bgn+v(w))$

and a node is a leaf if the condition of $n$ for sum is $a_{w,B}(n)=k$, has a
child corresponding to the right side if we are in the first case in Lemma 2,
and has $B$ children corresponding to the $B$ terms (minus signs are included
in the terms) of the right side if we are in the second case. Then $b_{w}(n)$
is the sum of the summands of the leaves of this tree. To prove (2.6), we
first notice that

$\log B\cdot
L(an+b)=\log\left(1-\frac{1}{an+b+1}\right)=-\frac{1}{an}+\mathcal{O}(1/n^{2})$

where $a$ and $b$ are positive constants. In particular, in (2.1),

$L(Bgn+v(w))=-\frac{1}{\log B\cdot B^{|w|}n}+\mathcal{O}(1/n^{2}).$

Then, we note that, in $T$, the first order term in the summand of each node
that is not a leaf is the sum of the first order term in its children. For
example, when a node has $B$ children and the first order term of the summand
is $-\dfrac{1}{\log B\cdot B^{s}n}$, then the sum of first order terms of the
summands of its children is

$-\frac{1}{\log B\cdot B^{s-1}n}-\sum_{\begin{subarray}{c}j=0\\\ j\neq
b_{1}\end{subarray}}^{B-1}\left(-\frac{1}{\log B\cdot
B^{s}n}\right)=-\frac{1}{\log B\cdot B^{s}n}.$

By induction we conclude that

$\log_{B}(b_{w}(n))=-\frac{1}{\log B\cdot B^{|w|}n}+\mathcal{O}(1/n^{2}).\qed$

######  Remark 1.

Theorem 3 above generalizes the case $B=2$ in [4] (also see [5]). It can also
give another proof of [19, Theorem 3].

######  Remark 2.

Actually the same “reducing trick” can be used to re-prove Equation (1.1) by
using a result in [3]. Namely, up to notation, it was proved in [3, Lemme, p.
142] that, for all $k\geq 0$,

$\sum_{s_{B}(n)=k}\log\left(\frac{n+1}{B\lfloor n/B\rfloor+B}\right)=-\log B.$
(2.7)

Define the fractional part of $n/B$ by $\\{n/B\\}:=n/B-\lfloor n/B\rfloor$.
Then, we have when $n$ tends to infinity,

$\begin{array}[]{lll}\displaystyle\log\left(\frac{n+1}{B\lfloor
n/B\rfloor+B}\right)&=&\displaystyle\log\left(1+\frac{1-B+B\\{n/B\\}}{n+B(1-\\{n/B\\})}\right)\\\
&=&\displaystyle\log\left(1+\frac{1-B+B\\{n/B\\}}{n}\right)+O(1/n^{2})\\\
&=&\displaystyle\frac{1-B}{n}+\left(\frac{B\\{n/B\\}}{n}\right)+O(1/n^{2}).\\\
\end{array}$

Thus (convergences are consequences of, e.g., Equation 2.7, see [3]):

$\displaystyle-\log B$
$\displaystyle=\sum_{s_{B}(n)=k}\log\left(\frac{n+1}{B\lfloor
n/B\rfloor+B}\right)$ $\displaystyle=(1-B)\sum_{s_{B}(n)=k}\frac{1}{n}\
+\sum_{s_{B}(n)=k}\frac{B\\{n/B\\}}{n}+O(1/n^{2}).$

Hence (note that the term $O(1/n^{2})$ below can be chosen independent of
$k$),

$\displaystyle-\log B$ $\displaystyle=(1-B)\sum_{s_{B}(n)=k}\frac{1}{n}+\
\sum_{0\leq j\leq B-1}\sum_{s_{B}(n)=k-j}\frac{j}{Bn+j}\ +O(1/n^{2})$
$\displaystyle=(1-B)\sum_{s_{B}(n)=k}\frac{1}{n}+\ \sum_{0\leq j\leq
B-1}\sum_{s_{B}(n)=k-j}\frac{j}{Bn}\ +O(1/n^{2}).$

Now, if $k$ tends to infinity, we have that $k-j$ tends to infinity for
$j\in[0,B-1]$, and also that $n$ must tend to infinity, hence, letting
$\lim_{k\to\infty}\sum_{s_{B}(n)=k}\frac{1}{n}:=\ell$,

$-\log B=(1-B)\ \ell+\sum_{0\leq j\leq B-1}\frac{j}{B}\ \ell=-\frac{B-1}{2}\
\ell,\ \ \ \text{thus \ \ $\ell=\frac{2\log B}{B-1}\cdot$}$

##  3 Statements and Declarations

The authors have no competing interests. No funds, grants, or other support
were received.

## References

  * [1]
  * [2] Alexander, R.: Remarks about the digits of integers. J. Austral. Math. Soc. 12, 239–241 (1971)
  * [3] Allouche, J.-P., Cohen, H., Mendès France, M., Shallit, J.: De nouveaux curieux produits infinis. Acta Arith. 49, 141–153 (1987)
  * [4] Allouche, J.-P., Shallit, J. O.: Infinite products associated with counting blocks in binary strings. J. London Math. Soc. 39, 193–204 (1989)
  * [5] Allouche, J.-P., Hajnal, P., Shallit, J. O.: Analysis of an infinite product algorithm. SIAM J. Discrete Math. 2, 1–15 (1989)
  * [6] Allouche, J.-P., Morin, C.: Kempner-like harmonic series. Preprint at https://arxiv.org/abs/2305.18180 (2023)
  * [7] Aloui, K.: Sur les entiers ellipséphiques : somme des chiffres et répartition dans les classes de congruence. Period. Math. Hungar. 70, 171–208 (2015)
  * [8] Aloui, K., Mauduit, C., Mkaouar, M.: Somme des chiffres et répartition dans les classes de congruence pour les palindromes ellipséphiques. Acta Math. Hungar. 151, 409–455 (2017)
  * [9] Baillie, R.: Sums of reciprocals of integers missing a given digit. Amer. Math. Monthly 86, 372–374 (1979). Also see ERRATA. Amer. Math. Monthly 87, 866 (1980)
  * [10] Baillie, R.: Summing the curious series of Kempner and Irwin. Preprint at https://arxiv.org/abs/0806.4410 (2023)
  * [11] Behforooz, G. H.: Thinning out the harmonic series. Math. Mag. 68, 289–293 (1995)
  * [12] Boas, R. P.: Some remarkable sequences of integers. In: Honsberger, R. (ed.) Mathematical Plums. The Dolciani Mathematical Expositions, 4, pp. 38–61. Mathematical Association of America, Washington, D.C. (1979)
  * [13] Col, S.: Propriétés multiplicatives d’entiers soumis à des conditions digitales. Thèse, Nancy (2006)
  * [14] Craven, B. D.: On digital distribution in some integer sequences. J. Austral. Math. Soc. 5, 325–330 (1965)
  * [15] Farhi, B.: A curious result related to Kempner’s series. Amer. Math. Monthly 115, 933–938 (2008)
  * [16] Fischer, H.-J.: Die Summe der Reziproken der natürlichen Zahlen ohne Ziffer $9$. Elem. Math. 48, 100–106 (1993)
  * [17] Gordon, R. A.: Comments on “Subsums of the harmonic series”. Amer. Math. Monthly 126, 275–279 (2019)
  * [18] Hu, N.: Fractal Uncertainty Principles for Ellipsephic Sets. M.Sc. Thesis, University of British Columbia (2021) http://dx.doi.org/10.14288/1.0396939
  * [19] Hu, Y.: Patterns in numbers and infinite sums and products. J. Number Theory 162, 589–600 (2016)
  * [20] Irwin, F.: A curious convergent series. Amer. Math. Monthly 23, 149–152 (1916)
  * [21] Kempner, A. J.: A curious convergent series. Amer. Math. Monthly 21, 48–50 (1914)
  * [22] Kløve, T.: Power sums of integers with missing digits. Math. Scand. 28, 247–251 (1971)
  * [23] Köhler, G., Spilker, J.: Dirichlet-Reihen zu Kempners merkwürdiger konvergenter Reihe. Math. Semesterber. 56, 187–199 (2009)
  * [24] Lubeck, B., Ponomarenko, V.: Subsums of the harmonic series. Amer. Math. Monthly 125, 351–355 (2018)
  * [25] Mukherjee, R., Sarkar, N.: A short note on a curious convergent series. Asian-Eur. J. Math. 14, Paper No. 2150158 (2021)
  * [26] Nathanson, M. B.: Dirichlet series of integers with missing digits. J. Number Theory 222, 30–37 (2021)
  * [27] Nathanson, M. B.: Convergent series of integers with missing digits. Ramanujan J. 58, 667–676 (2022)
  * [28] Nathanson, M. B.: Curious convergent series of integers with missing digits. Integers 21A, Ron Graham Memorial Volume, Paper No. A18 (2021)
  * [29] Schmelzer, T., Baillie, R.: Summing a curious, slowly convergent series. Amer. Math. Monthly 115, 525–540 (2008)
  * [30] Segal, A. C., Lepp, B., Fine, N. J.: A limit problem [E 2204]. Amer. Math. Monthly 77, 1009–1010 (1970)
  * [31] Wadhwa, A. D.: An interesting subseries of the harmonic series. Amer. Math. Monthly 82, 931–933 (1975)
  * [32] Wadhwa, A. D.: Some convergent subseries of the harmonic series. Amer. Math. Monthly 85, 661–663 (1978)
  * [33] A. Walker, A. Walker, Arithmetic progressions with restricted digits. Amer. Math. Monthly 127, 140–150 (2020)


# A preference elicitation interface for collecting dense recommender datasets
with rich user information

Demo

Pantelis P. Analytis  Cornell University<EMAIL_ADDRESS>,  Tobias Schnabel
Cornell University<EMAIL_ADDRESS>,  Stefan Herzog  MPI for Human
Development<EMAIL_ADDRESS>,  Daniel Barkoczi  MPI for Human
Development<EMAIL_ADDRESS>and  Thorsten Joachims  Cornell
University<EMAIL_ADDRESS>
###### Abstract.

We present an interface that can be leveraged to quickly and effortlessly
elicit people’s preferences for visual stimuli, such as photographs, visual
art and screensavers, along with rich side-information about its users. We
plan to employ the new interface to collect dense recommender datasets that
will complement existing sparse industry-scale datasets. The new interface and
the collected datasets are intended to foster integration of research in
recommender systems with research in social and behavioral sciences. For
instance, we will use the datasets to assess the diversity of human
preferences in different domains of visual experience. Further, using the
datasets we will be able to measure crucial psychological effects, such as
preference consistency, scale acuity and anchoring biases. Last, we the
datasets will facilitate evaluation in counterfactual learning experiments.

preference elicitation, recommender system datasets, visual art

††copyright: none††ccs: Human-centered computing Collaborative filtering††ccs:
Human-centered computing Social media††ccs: Human-centered computing
Collaborative and social computing devices

##  1\. Introduction

Over the last three decades the recommender systems community has made immense
progress in the way we represent, understand and learn people’s preferences as
a function of previously collected explicit or implicit evaluations. Research
in recommender systems has by all means increased the quality of the curated
and recommended content in the online world. Several large datasets have been
a crucial component of this success, as they have commonly functioned as test-
beds on which new theories and algorithms have been compared (Movielens,
LastFM and Netflix to name just a few). Most of these datasets, however, are
very sparse. They contain thousands items and even the most popular among the
items have been evaluated only by a small subset of their users. Given the
large fraction of missing ratings, it is challenging to accurately estimate
even simple quantities like the average quality of an item, especially since
the patterns of missing data are subject to strong selection biases
(Pradel/etal/12, ). This presents fundamental challenges when evaluating
recommendation algorithms on sparse datasets. Further, it becomes an obstacle
for scholars in the social and behavioral sciences as workarounds have to be
developed for dealing with missing values.

To the best of our knowledge, the only dense collaborative filtering dataset
was the outcome of the Jester Interface (goldberg2001eigentaste, ). The
interface curated 100 jokes of various styles and topics. People utilized a
slider to evaluate 5 jokes that were presented to them sequentially. The first
evaluations were used to estimate people’s preferences and to recommend them
the remaining jokes. The users continued to read and evaluate jokes until the
pool of 100 items was exhausted. In total, more than 70.000 people have
evaluated at least some of the jokes, and more than 14.000 have evaluated all
the jokes, resulting to a fully evaluated subset of the dataset.

Figure 1. The design of the preference elicitation interface. We replicate the
design of the _Jester_ interface, using a continuous bar that people can use
to express how much they liked or disliked an item. Participants have to wait
for at least 5 seconds before they can proceed to the next item.

##  2\. The interface and data collection

We plan to collect new datasets in different domains of people’s visual
experience, ranging from photographs and paintings to designs for
screensavers. Our interface replicates the design of the Jester interface,
adding new elements that can counteract its limitations. At the outset, people
are provided with instructions about how to use the interface. Then, before
the presentation of the stimuli we collect demographic information about the
users. To reduce possible order effects, the visual stimuli are presented in
random order. As in Jester, users are asked to evaluate items using a slider
bar; they can move the marker of the slide bar to the left to indicate that
they did not like the item, or to the right to indicate that they liked it. We
implement a continuous scale, which allows a fine-grained evaluation of the
presented items. Finally, to limit anchoring bias, the slide bar is initially
semi-transparent and the colors become vivid only when the user has clicked on
it.111The interface can be accessed at http://abc-webstudy.mpib-
berlin.mpg.de/recstrgs/study_simulator.php. Both the code for the interface
and the collected data will be publicly available.

Once all the items have been evaluated, we collect further psychologically
relevant information about the users. Numerous studies have shown that side
information can substantially improve estimates of people’s preference and it
complements first hand evaluations (park2009pairwise, ). In the first
experiments we will deploy the visual-art expertise questionnaire developed by
Chatterjee et al. (chatterjee2010assessment, ) to gauge people’s familiarity
with the visual arts and a succinct version of the big-five questionnaire to
quickly assess the people’s personalities (rammstedt2007measuring, ) (see
Figure 2). It takes about 20 minutes to complete the current version of the
interface, including the instructions, questionnaires and evaluation phase.

We intend to conduct the first experiments at Amazon’s Mechanical Turk labor
market. Several studies have shown that for effortless tasks the results
produced on mTurk are comparable to laboratory studies (paolacci2010running,
). The visual stimuli used in this interface evoke immediate aesthetic
judgments, and thus can quickly be transformed to evaluations. Eventually, we
intend to develop a data visualization tool that will reward people who
complete the study with information about their preference profiles and how
they relate to those of other individuals. Thus, we intend to create an
inherently motivating interface using as a reward the informational value
generated by the collected data. In this way, we will reduce the cost of data
collection, but also introduce basic ideas behind collaborative filtering and
recommender systems to the wider public.

Figure 2. At the end of the evaluation phase we collect additional information
about the users. We invited the users to complete a questionnaire about their
expertist in the visual arts and a 10-question version of the big-five
questionnaire.

##  3\. Potential applications

We envisage several new applications for the developed datasets. Here we
foreshadow a few of these potential uses, keeping in mind that the community
that will have access to the produced datasets will certainly come up with
more. First, they will facilitate cross-fertilization with the cognitive and
behavioral sciences. For instance, social and cognitive psychologists have
extensively studied simple strategies for inference and estimation where
different features are used to predict an objective truth. The new datasets
will open the way to study strategies for social preference learning in
domains where no objective truth exists (analytis2015you, ). Also, we can
manipulate the design of the interface to study relevant behavioral effects,
such as to study the consistency of evaluations or to investigate the effect
of the granularity of the evaluation scale on the predictions. To sum up, the
datasets will allow us to better understand preference diversity and its
implications for different recommender systems algorithms as well as for
psychological social learning strategies.

Moreover, we believe that the new datasets can fuel existing streams of
research in recommender systems and machine learning. For instance, dealing
with selection-biases and with data missing not at random is a growing
research stream in recommender systems and machine learning
(schnabel2016unbiased, ). To evaluate algorithms tuned to deal with such
problems, we can impose selection biases ex-ante and remove data from the
dense dataset accordingly. This set-up could complement existing sparse
datasets for learning, with the difference that selection biases can be
controlled and varied in order to test robustness. Moving on to the broader
class of counterfactual simulations, dense datasets greatly simplify
evaluation since they can serve as ground-truth when conducting simulations
(salganik2006experimental, ).

## References

  * [1] Pantelis P Analytis, Daniel Barkoczi, and Stefan M Herzog.  You’re special, but it doesn’t matter if you’re a greenhorn: Social recommender strategies for mere mortals.  In Cognitive Science Society, pages 1799–1804. Cognitive Science Society, 2015.
  * [2] Anjan Chatterjee, Page Widick, Rebecca Sternschein, William B Smith, and Bianca Bromberger.  The assessment of art attributes.  Empirical Studies of the Arts, 28(2):207–222, 2010.
  * [3] Ken Goldberg, Theresa Roeder, Dhruv Gupta, and Chris Perkins.  Eigentaste: A constant time collaborative filtering algorithm.  Information Retrieval, 4(2):133–151, 2001.
  * [4] Gabriele Paolacci, Jesse Chandler, and Panagiotis G Ipeirotis.  Running experiments on Amazon Mechanical Turk.  Judgment and Decision Making, 5(5), 2010.
  * [5] Seung-Taek Park and Wei Chu.  Pairwise preference regression for cold-start recommendation.  In Proceedings of the third ACM conference on Recommender systems, pages 21–28. ACM, 2009.
  * [6] B. Pradel, N. Usunier, and P. Gallinari.  Ranking with non-random missing ratings: influence of popularity and positivity on evaluation metrics.  In RecSys, pages 147–154, 2012.
  * [7] Beatrice Rammstedt and Oliver P John.  Measuring personality in one minute or less: A 10-item short version of the big five inventory in english and german.  Journal of Research in Personality, 41(1):203–212, 2007.
  * [8] Matthew J Salganik, Peter Sheridan Dodds, and Duncan J Watts.  Experimental study of inequality and unpredictability in an artificial cultural market.  Science, 311(5762):854–856, 2006.
  * [9] Tobias Schnabel, Adith Swaminathan, Peter I Frazier, and Thorsten Joachims.  Unbiased comparative evaluation of ranking functions.  In ICTIR, pages 109–118, 2016.


# An analytical reconstruction formula with efficient implementation for a
modality of Compton Scattering Tomography with translational geometry

Cécilia Tarpaua,b,c,∗ ID , Javier Cebeirod ID , Geneviève Rolleta ID ,
Mai K. Nguyenb ID and Laurent Dumasc ID
a LPTM (UMR 8089), CY Cergy Paris Université, CNRS, Cergy-Pontoise, France
bETIS (UMR 8051), CY Cergy Paris Université, ENSEA, CNRS, Cergy-Pontoise,
France
c LMV (UMR 8100), Université de Versailles Saint-Quentin, CNRS, Versailles,
France
d CEDEMA, Universidad Nacional de San Martín, San Martín, Argentina
∗ Corresponding author

###### Abstract

In this paper, we address an alternative formulation for the exact inverse
formula of the Radon transform on circle arcs arising in a modality of Compton
Scattering Tomography in translational geometry proposed by Webber and Miller
(Inverse Problems (36)2, 025007, 2020). The original study proposes a first
method of reconstruction, using the theory of Volterra integral equations. The
numerical realization of such a type of inverse formula may exhibit some
difficulties, mainly due to stability issues. Here, we provide a suitable
formulation for exact inversion that can be straightforwardly implemented in
the Fourier domain. Simulations are carried out to illustrate the efficiency
of the proposed reconstruction algorithm.

Keywords: Analytic inversion, Compton Scattering Tomography, Image formation,
Image reconstruction, Radon transform on double circle arcs

##  1 Introduction

Compton Scattering Tomography (CST) is an imaging technique whose objective is
to exploit wisely Compton scattered photons by the object to scan in order to
reconstruct its electron density map. Since the early proposition of this type
of imaging by Lale [1], Clarke [2] and Farmer [3], studies about CST systems
proved already promising results in medical imaging, for instance for the
identification of lung tumours [4, 5], but also for earthquake engineering
[6], cultural heritage imaging [7, 8], landmine detection [9] and agricultural
measurements [10]. In fact, CST has made it possible to widen the fields of
application for tomography to the imaging of one-sided large objects, because,
in some configurations, sources and detectors can be placed at the same side
of the object.

The proposition of CST systems is strongly related to the study of the
associated integral transform, which models data acquisition [11, 12]. These
integral transforms are generalizations of the classical Radon transform on
lines studied by Radon [13] and Cormack [14]. While, in two-dimensional CST,
manifolds are on circle arcs or double circle arcs, in three dimensions, data
is acquired on toric surfaces. These circular geometries for the considered
manifolds originate from the Compton effect. This physical phenomenon occurs
when a photon emitted by the source with energy $E_{0}$ collides with an
electron as it passes through matter. This photon is scattered, and deviated
by an angle $\omega$ from its original direction. The photon looses also a
part of its energy. The Compton formula gives us the one-to-one correspondence
between the energy of scattered photons $E(\omega)$ and the related scattering
angle $\omega$

$E(\omega)=\frac{E_{0}}{1+\frac{E_{0}}{m\,c^{2}}\left(1-\cos(\omega)\right)}.$
(1.1)

$m$ is the electron mass and $c$ the speed of light. This relation ensures
also that scattering sites of photons with identical energy $E(\omega)$ are
located on a circle arc labelled by the scattering angle $\omega$. Figure 1
illustrates the general functioning principle of a CST system.

Figure 1: General functioning principle of a CST system. Photons are emitted
by source $S$, interact at sites $M$, and are recorded at site $D$. When a
photon is detected carrying an energy $E(\omega_{1})$ (resp. $E(\omega_{2})$),
the possible interaction sites lie on the upper (resp. lower) circle arc which
subtends the angle $(\pi-\omega_{1})$ (resp. $(\pi-\omega_{2})$).

Several issues around the study of these generalized Radon transforms are then
of interest, such as existence and uniqueness of the solution, its stability
or the range conditions. The main important problem remains the reconstruction
of the image, mostly coming from the proposition of analytical inverse
formulas.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: Previous proposed CST modalities (a, b, c, d, e, f). (a): Fixed
source and detectors placed on a line. (b): Rotating pair source-detector
diametrically opposed. (c): Rotating pair source detector. (d): Fixed source
and detectors placed on a ring. (e): Detector rotating around a fixed source.
(f) Source-detector translating simultaneously along two parallel lines. In
all figures: The source $S$ is represented by a red point. The detector(s) $D$
is (are) represented by blue point(s). The $M$, $M_{i}$ or $M^{\prime}_{i}$ in
black, are running points and examples of scattering site. An example of
trajectory for a photon whose scattering site is $M$ is shown in purple. The
corresponding scattering angle is denoted $\omega$. The object to scan is
represented in grey. The red continuous curves are the examples of scanning
circles arcs. For (c), (e) and (f), the dashed circles (resp. lines)
represents the circular (resp. linear) paths on which move the sensors.

In that way, the first studied two-dimensional CST system was made of a fixed
line of detectors containing a fixed source [15]. See Figure 2a. The
associated scanning circle arcs consist of a family of semicircles with a
common fixed end point (the point source) and another extremity (the
considered detector) on the line. Then, circular geometries for CST have been
proposed, the first one considers a pair source - detector diametrically
opposed in rotation around the object [16, 17, 18] and data acquisition is
modelled by a Radon transform on circle arcs having a fixed source-detector
chord length. See Figure 2b. The second consists also of a pair source -
detector moving on a circle, but the distance between the two is no longer
constant and depends on the scattering angle [19] (Fig. 2c). With this
modality, data measurement is performed on a family of circles orthogonal to
the fixed circular path of the pair source - detector. In a third modality
(see Fig. 2d), on the contrary, it was considered a set of detectors placed on
a ring containing the source, thus obtaining a completely fixed CST modality
[20, 21, 22, 23]. The corresponding Radon transform measures the contribution
of the photons scattered at the points located on circles arcs having a common
extremity (the point source) and another (the considered detector) on the
detector ring. Some configurations, mentioned above, employ collimators to
split up photons coming from different circle arcs. Geometries without
collimation have also been studied. Consequently, for a given position of
detector, the acquisition is performed on a family of double circle arcs and
the amount of registered data is potentially doubled. This feature may also
reduce the acquisition time and finally radiation exposition in comparison
with a similar geometry with collimation. Among the proposed CST systems
without collimation, one supposes a fixed source and a detector rotating
around this source [24]. See Figure 2e. The other one, proposed in [25], is
made of a source and a detector which translate along a line (Fig. 2f).

The direct reconstruction of volumes is also of interest with the proposition
of three-dimensional CST systems. These modalities use uncollimated sources
and detectors, and data acquisition is performed on toric surfaces. In many
cases, these 3D modalities correspond to extensions of 2D systems. Circular
geometries become thus spherical or cylindrical [26, 27, 28] and linear
geometries become planar [25].

Here, we are interested in the two-dimensional modality proposed in [25]. The
purpose of this article is to develop an alternative formulation suitable for
the development of a faster and efficient reconstruction algorithm. The
associated reconstruction algorithm will use only classical tools such as Fast
Fourier Transform (FFT) algorithm.

The paper is outlined as follows. Section 2 recalls the general setup of the
system and the model for data acquisition for the modality. Section 3
introduces the main result of the paper, that is, the alternative formulation
for the inverse Radon transform on double circle arcs. Section 4 will give the
discrete formulation of the forward operator, as well as the proposed strategy
to reconstruct the object under study. Section 5 discusses the obtained
simulation results with a study of the influence of some parameters on
reconstruction quality.

##  2 Setup and measurement model of the CST system under study

###  2.1 Setup

The system under study is made of a source, assumed to be monochromatic, and a
detector separated by a fixed distance from each other. The source and the
detector move respectively on a horizontal line of equation $z=3$ and $z=1$.
The horizontal position of the pair source-detector is labelled by $x_{0}$
(see Figure 3). Alternatively, this system may be sketched with fixed lines of
sources and detectors that will be used in pair. The object, placed below the
detector path, is scanned transversely. No collimation is used at the
detector, hence the acquisition is performed on a family of double circle arcs
(called toric sections in the original publication [25])111We position
ourselves in the same frame of study as the original article, with the same
working assumptions. Thus, first order Compton scattering is the only source
of attenuation for radiation and data acquisition is performed with a pair
source detector, assumed to be point-like. These conditions are common in the
literature [16, 26, 29, 30, 31, 20, 24] and have been already discussed in
[15, 19, 26, 24].. We parameterize these circle arcs and define the
corresponding Radon transform in the next paragraph.

Figure 3: Setup and parameterization of the CST modality proposed in [25]. The
source $S$ and the detector are respectively represented by a red and a blue
point. To make the difference between the four half-arcs, $S_{1},S_{3}$ and
$S_{2},S_{4}$ are respectively depicted in red and green. $\Omega_{1,2,3,4}$
denote the centres of the circles supporting the half-arcs $S_{1,2,3,4}$. The
point $M$ is an example of a scattering site.

###  2.2 Modelling of data acquisition using the CST system

Given a scattering angle $\omega$, data acquisition is performed on a family
of double circle arcs of radius $r$ where $r=1/\sin(\pi-\omega)$ (or,
equivalently, $\omega=\pi-\arcsin{(1/r)}$). For parameterization, these double
circle arcs are obtained with the union of four half arcs denoted
$S_{j}(x_{0},r)$, $j\in\\{1,2,3,4\\}$ of respective equation

$\displaystyle x_{1}$
$\displaystyle=x_{0}+\sqrt{r^{2}-1}+\sqrt{r^{2}-(z-2)^{2}},\;\;x_{2}=x_{0}+\sqrt{r^{2}-1}-\sqrt{r^{2}-(z-2)^{2}},$
$\displaystyle x_{3}$
$\displaystyle=x_{0}-\sqrt{r^{2}-1}+\sqrt{r^{2}-(z-2)^{2}},\;\;x_{4}=x_{0}-\sqrt{r^{2}-1}-\sqrt{r^{2}-(z-2)^{2}}$

and $z\in]2-r,1[$. See Figure 3.

The Radon transform which mathematically models data measurement with this CST
system is then defined as follows :

######  Definition 1.

Let $f$ be an unknown function, non-negative, continuous and compactly
supported in the half plane $z<1$. The Radon transform on double circle arcs
$\mathcal{R}_{\mathcal{D}}$ maps $f$ into the set of its integrals over the
family of double circle arcs as

$\mathcal{R}_{\mathcal{D}}f(x_{0},r)=\int_{\bigcup_{j=1}^{4}S_{j}(x_{0},r)}f(x,y)ds.$
(2.1)

where $ds$ refers to the elementary arc length measure on the considered
double circle arc. Then, after computation of the arc length measure, we have
the explicit reformulation for $\mathcal{R}_{\mathcal{D}}$ [25, Proposition
3.1] :

$\mathcal{R}_{\mathcal{D}}f(x_{0},r)=\int_{1}^{r}\frac{1}{\sqrt{1-\left(\frac{z}{r}\right)^{2}}}\left(\sum_{j=1}^{2}f_{1}\left(\sqrt{r^{2}-1}+(-1)^{j}r\sqrt{1-\left(\frac{z}{r}\right)^{2}}+x_{0},z\right)\right.+\\\
\left.f_{1}\left(-\sqrt{r^{2}-1}+(-1)^{j}r\sqrt{1-\left(\frac{z}{r}\right)^{2}}+x_{0},z\right)\right)dz,$
(2.2)

where $f_{1}(x,z)=f(x,2-z)$.

In the original study of this modality, the invertibility of the corresponding
Radon transform as well as its analytical inversion formula has been
established. The invertibility was proven using the theory of integral
equations and resulted in a Volterra integral equation with a weakly singular
kernel in the Fourier domain. This study also leads to a formulation for
inversion formula as an integral transformation with a kernel computed
iteratively. The numerical calculation of this kind of kernel may require high
computational time and/or memory. Furthermore, as mentioned in the Remark 3.4
of the original paper, the proposed approach by Webber is severely ill-posed,
particularly in terms of stability. Implementing such a method can lead to
large instabilities, even when these are due to small changes in the data.

##  3 An alternative formulation for the inversion formula of the Radon
transform on double circle arcs

In this section, we state the main result of the paper, a different
formulation for the associated inversion formula, that will be easier to
implement numerically. Let us introduce before some notations that will be
used in the proofs.

###  3.1 Notations

It is useful to define the following transform pairs.

######  Definition 2 (Fourier transform).

Let $f$ be a compactly supported function in $\mathbb{R}^{n}$. The
n-dimensional Fourier transform of $f$, denoted $\widehat{f}$, is given by

$\widehat{f}(\boldsymbol{\xi})=\int_{\mathbb{R}^{n}}f(x)e^{-i\boldsymbol{x}\cdot\boldsymbol{\xi}}d\boldsymbol{x}$
(3.1)

with $\boldsymbol{\xi}\in\mathbb{R}^{n}$. The inverse Fourier transform is

$f(\boldsymbol{x})=\frac{1}{(2\pi)^{n}}\int_{\mathbb{R}^{n}}\widehat{f}(\boldsymbol{\xi})e^{i\boldsymbol{x}\cdot\boldsymbol{\xi}}d\boldsymbol{\xi}.$
(3.2)

######  Definition 3 (Fourier cosine transform [32]).

Let $f$ be a compactly supported function in $\mathbb{R}^{+}$. The Fourier
cosine transform of $f$, denoted $\widehat{f}^{c}$, is given by

$\widehat{f}^{c}(\xi)=\sqrt{\frac{2}{\pi}}\int_{0}^{\infty}f(x)\cos(x\xi)dx$
(3.3)

with $\xi\in\mathbb{R}$. The inverse Fourier cosine transform is

$f(\xi)=\sqrt{\frac{2}{\pi}}\int_{0}^{\infty}\widehat{f}^{c}(\xi)\cos(x\xi)d\xi.$
(3.4)

We define also the Hankel transform.

######  Definition 4 (Hankel transform [33]).

Let $f$ be a compactly supported function in $\mathbb{R}^{+}$. The zero-order
Hankel transform of $f$ is defined as

$\mathcal{H}_{0}f(\eta)=\int_{0}^{\infty}f(r)J_{0}(\eta r)rdr$ (3.5)

where $J_{0}$ stands for the Bessel function of the first kind of order $0$.

Finally, we recall the integral representation of the Bessel function $J_{0}$:

$J_{0}(x)=\frac{1}{2\pi}\int_{-\pi}^{\pi}e^{ix\sin{\theta}}d\theta.$ (3.6)

###  3.2 Inversion formula

######  Proposition 1.

Denoting $\mathcal{G}f(x_{0},r)$ the operator whose Fourier transform
according to the first variable is

$\widehat{\mathcal{G}}f(\xi,r)=\frac{\widehat{\mathcal{R}}_{\mathcal{D}}f(\xi,r)}{2r\cos{(\xi\sqrt{r^{2}-1})}},$
(3.7)

if $r>1$ and $0$ when $r\in[0,1]$, the unknown function $f$ is completely
recovered from $\widehat{\mathcal{G}}f$ as follows

$f(x,z)=\frac{1}{4\pi}\int_{-\infty}^{\infty}e^{ix\xi}\int_{0}^{\infty}\mathcal{H}_{0}\widehat{\mathcal{G}}f(\xi,\sqrt{\xi^{2}+\sigma^{2}})\cos{(\sigma(2-z))}\sigma
d\sigma d\xi.$ (3.8)

###### Proof.

With the change of variables $s=z/r$ in (2.2) and taking the Fourier transform
of $\mathcal{R}_{\mathcal{D}}$ respectively to variable $x_{0}$, one gets

$\widehat{\mathcal{R}}_{\mathcal{D}}f(\xi,r)=\int_{-\infty}^{\infty}dx_{0}\int_{1/r}^{1}ds\frac{re^{-ix_{0}\xi}}{\sqrt{1-s^{2}}}\cdot\\\
\left(\sum_{j=1}^{2}f_{1}\left(\sqrt{r^{2}-1}+(-1)^{j}r\sqrt{1-s^{2}}+x_{0},rs\right)+f_{1}\left(-\sqrt{r^{2}-1}+(-1)^{j}r\sqrt{1-s^{2}}+x_{0},rs\right)\right).$
(3.9)

With the second change of variables
$x=x_{0}\pm\sqrt{r^{2}-1}+(-1)^{j}r\sqrt{1-s^{2}}$, one gets

$\displaystyle\widehat{\mathcal{R}}_{\mathcal{D}}f(\xi,r)=4\int_{1/r}^{1}ds\frac{r}{\sqrt{1-s^{2}}}\widehat{f}_{1}(\xi,rs)\cos{(\xi\sqrt{r^{2}-1})}\cos{(\xi
r\sqrt{1-s^{2}})}$ (3.10)

where $\widehat{f}_{1}$ stands for the one-dimensional Fourier transform
relatively to the variable $x$.

Using relation (3.7), multiplying both sides of (3.10) by $r\cdot J_{0}(\eta
r)$ with $\eta\geq 1$ and integrating with respect to variable $r$, for $r>1$,
one recognizes the Hankel transform of $\widehat{\mathcal{G}f}$, denoted
$\mathcal{H}_{0}\widehat{\mathcal{G}f}$

$\mathcal{H}_{0}\widehat{\mathcal{G}}f(\xi,\eta)=2\int_{1}^{\infty}dr\int_{1/r}^{1}ds\frac{r}{\sqrt{1-s^{2}}}\widehat{f}_{1}(\xi,rs)\cos{(\xi
r\sqrt{1-s^{2}})}J_{0}(\eta r).$ (3.11)

Then, with the double substitution
$(r=\sqrt{z^{2}+b^{2}},s=z/\sqrt{z^{2}+b^{2}})$, one gets

$\mathcal{H}_{0}\widehat{\mathcal{G}}f(\xi,\eta)=2\int_{1}^{\infty}dz\widehat{f}_{1}(\xi,z)\int_{0}^{\infty}db\cos{(\xi
b)}J_{0}(\eta\sqrt{z^{2}+b^{2}}).$ (3.12)

The result of the $b$-integral is given in the table [34, p. 55, eq. (35)].
Finally, one gets for $0<\xi<\eta$

$\mathcal{H}_{0}\widehat{\mathcal{G}}f(\xi,\eta)=2\int_{1}^{\infty}dz\widehat{f}_{1}(\xi,z)\frac{1}{\sqrt{\eta^{2}-\xi^{2}}}\cos{(z\sqrt{\eta^{2}-\xi^{2}})}$
(3.13)

and $\mathcal{H}_{0}\widehat{\mathcal{G}}f(\xi,\eta)=0$ if $\eta\leq\xi$.

Let $0<\xi<\eta$. Then with the fact $f_{1}(x,z)=0$ for $z\in[0,1]$,

$\int_{0}^{\infty}dz\widehat{f}_{1}(\xi,z)\cos{(z\sqrt{\eta^{2}-\xi^{2}})}=\frac{\sqrt{\eta^{2}-\xi^{2}}}{2}\mathcal{H}_{0}\widehat{\mathcal{G}}f(\xi,\eta).$
(3.14)

The left-hand side is the Fourier cosine transform of $\widehat{f}_{1}(\xi,z)$
according the variable $z$. We can then extract $\widehat{f}_{1}(\xi,z)$,
applying the inverse cosine transform to (3.14)

$\widehat{f}_{1}(\xi,z)=\frac{1}{2}\int_{0}^{\infty}d\sigma\mathcal{H}_{0}\widehat{\mathcal{G}}f(\xi,\sqrt{\xi^{2}+\sigma^{2}})\cos{(z\sigma)}\sigma$
(3.15)

where $\sigma=\sqrt{\eta^{2}-\xi^{2}}$. The final equation is obtained going
back to variable $f$, and applying the inverse Fourier transform. ∎

######  Remark 1.

The projections $\mathcal{\widehat{G}}f$ (3.7) contain zeros in the
denominator, since the cosine function vanishes when
$\xi\sqrt{r^{2}-1}=2k\pi\pm\frac{\pi}{2}$, $k\in\mathbb{Z}$. From (3.11) to
the end of the demonstration, it was supposed that $r$ is different from
$\sqrt{1+\left(\frac{\pi}{2\xi}(2k+1)\right)}$. Furthermore, this may be a
source of instability in the simulations. The addition of a regularization
parameter for simulations is discussed in Section 4.2 to prevent this.

######  Remark 2.

Another reconstruction algorithm is also possible from the projections of
$\mathcal{G}f$. This process is achieved performing geometric inversion.
Geometric inversion is a mapping converting a point $X$ into a point
$\tilde{X}$ such that $\tilde{X}X^{T}=q^{2}$, where $q\in\mathbb{R}_{+}^{*}$
is a constant value. The mapped point $\tilde{X}$ has the same direction as
the original point $X$ but a distance of $q^{2}/||X||$ to the origin of the
considered coordinate system. As an example, geometric inversion converts
circles passing through the origin into straight lines. In the present case,
the Radon transform on double circle arcs is converted into a Radon transform
on an apparent family of circle arcs of similar geometry as the one studied in
[19, 35]. Although the inverse problem can be alternatively solved using
geometric inversion, the approach we employ here is more straightforward.

##  4 Numerical formulations for the forward and inverse transform

###  4.1 Image formation

Let $N_{SD}$ be the number of positions for the pair source - detector and
$N_{r}$ the number of double scanning circle arcs per sensor position. We
denote $x_{0,k}$, $k\in\\{1,...,N_{SD}\\}$ and $r_{l}$,
$l\in\\{1,...,N_{r}\\}$ the discrete variables corresponding respectively to
$x_{0}$ and $r$. The matrix of projection data
$\mathcal{R}_{\mathcal{D}}f(x_{0,k},r_{l})$ is then computed, writing (2.2)
under a discrete form, with the change of variables $z=r\cos{\theta}$

$\mathcal{R}_{\mathcal{D}}f(x_{0,k},r_{l})=r_{l}\,\Delta_{\theta}\;\cdot\\\
\sum_{\theta\in\left[\arcsin{\left(\frac{1}{r_{l}}\right)},\frac{\pi}{2}\right]}\left(\sum_{j=1}^{2}f_{1}(x_{0,k}+\sqrt{r_{l}^{2}-1}+(-1)^{j}r_{l}\cos{\theta},r_{l}\sin{\theta})+f_{1}(x_{0,k}-\sqrt{r_{l}^{2}-1}-(-1)^{j}r_{l}\cos{\theta},r_{l}\sin{\theta})\right),$
(4.1)

where $\Delta_{\theta}$ is the sampling angular distance of $\theta$. The
above Cartesian parameterization allows having a constant distance between
running points of the considered scanning circle arcs during simulations.

###  4.2 Image reconstruction

For image reconstruction, we need to compute the projections $\mathcal{G}f$ in
the Fourier domain according to (3.7). This expression contains zeros in the
denominator. This may induce instabilities on reconstruction. For simulations,
we add a small regularization parameter denoted $\epsilon$

$\widehat{\mathcal{G}}f(\xi,r)=\frac{\widehat{\mathcal{R}}_{\mathcal{D}}f(\xi,r)}{2r}\frac{\cos{(\xi\sqrt{r^{2}-1})}}{\epsilon^{2}+\cos{(\xi\sqrt{r^{2}-1})}^{2}}.$
(4.2)

In terms of computational cost, the most demanding step in the implementation
of (3.8) is the calculation of $\mathcal{H}_{0}\widehat{\mathcal{G}}f$. The
idea is to establish a relation between the above operator with the Fourier
transform of $\mathcal{G}^{\ddagger}f$, defined as follows

$\mathcal{G}^{\ddagger}f(x,z)=\int_{-\infty}^{\infty}\mathcal{G}f(x_{0},\sqrt{(x-x_{0})^{2}+z^{2}})\,dx_{0}.$
(4.3)

We have now the following proposition.

######  Proposition 2.

Let $(\sigma,\xi)\in[0,\infty[\times\mathbb{R}$.
$\mathcal{H}_{0}\widehat{\mathcal{G}}f(\xi,\sqrt{\xi^{2}+\sigma^{2}})$ is
related with the two-dimensional Fourier transform of the operator
$\mathcal{G}^{\ddagger}$ as

$\mathcal{H}_{0}\widehat{\mathcal{G}}f(\xi,\sqrt{\xi^{2}+\sigma^{2}})=2\pi\widehat{\mathcal{G}}^{\ddagger}f(\xi,\sigma).$
(4.4)

Consequently, from the inversion formula (3.8), it follows in the Fourier
domain

$\widehat{f}_{1}(\xi,\sigma)=2\pi^{2}|\sigma|\widehat{\mathcal{G}}^{\ddagger}f(\xi,\sigma).$
(4.5)

where $\widehat{f}_{1}$ is the two-dimensional Fourier transform of $f_{1}$.

###### Proof.

From the definitions of Fourier and Hankel transforms and with the integral
representation of the Bessel function, one gets

$\displaystyle\mathcal{H}_{0}\widehat{\mathcal{G}}f(\xi,\sqrt{\xi^{2}+\sigma^{2}})$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}dx_{0}\int_{0}^{\infty}dr\,\mathcal{G}f(x_{0},r)\left(\int_{-\pi}^{\pi}e^{ir\sqrt{\xi^{2}+\sigma^{2}}\sin(\theta)}d\theta\right)\,e^{-i\xi
x_{0}}.$ (4.6)

There is an angle $\phi{(\sigma,\xi)}\in[0,2\pi[$ which corresponds to the
angular coordinate of the point $(\sigma,\xi)$, such that
$\sigma=\sqrt{\xi^{2}+\sigma^{2}}\cos{(\phi{(\sigma,\xi)})}$ and
$\xi=\sqrt{\xi^{2}+\sigma^{2}}\sin{(\phi{(\sigma,\xi)})}$. Using the property
of periodicity of trigonometric functions, it follows that

$\mathcal{H}_{0}\widehat{\mathcal{G}}f(\xi,\sqrt{\xi^{2}+\sigma^{2}})=\frac{1}{2\pi}\int_{-\infty}^{\infty}dx_{0}\int_{0}^{\infty}dr\,\mathcal{G}f(x_{0},r)\left(\int_{-\pi}^{\pi}e^{ir\xi\sin{\theta}+\sigma\cos{\theta}}d\theta\right)\,e^{-i\xi
x_{0}}.$ (4.7)

Changing variables $x=r\cos{\theta}$ and $z=r\sin{\theta}$,

$\mathcal{H}_{0}\widehat{\mathcal{G}}f(\xi,\sqrt{\xi^{2}+\sigma^{2}})=\frac{1}{2\pi}\int_{\mathbb{R}^{2}}dxdz\left(\int_{-\infty}^{\infty}dx_{0}\mathcal{G}f(x_{0},\sqrt{(x-x_{0})^{2}+z^{2}})\right)e^{-i(x\xi+z\sigma)}$
(4.8)

The right-hand side of (4.8) is the two-dimensional Fourier transform of
$\mathcal{G}^{\ddagger}$, weighted by $2\pi$ (4.4). We are now able to
reformulate the inversion formula (3.15) as

$\widehat{f}_{1}(\xi,z)=2\pi\int_{-\infty}^{\infty}d\sigma\widehat{\mathcal{G}}^{\ddagger}f(\xi,\sigma)e^{iz\sigma}|\sigma|.$
(4.9)

Taking the Fourier transform according to variable $z$ to the above equation
(4.9) leads to (4.5). ∎

This leads to the reconstruction algorithm summed up in Algorithm 1.

Data: $\mathcal{R}_{\mathcal{D}}f(x_{0},r)$, projections on double circular
arcs of function $f$

Result: $f(x,y)$

1 Compute the one-dimensional Fourier transform of
$\mathcal{R}_{\mathcal{D}}f(x_{0},r)$ using FFT relative to the first variable
;

2 Compute $\widehat{\mathcal{G}}f(\xi,r)$ according to (4.2) and perform the
inverse FFT to recover $\mathcal{G}f(x_{0},r)$;

3 For each $x_{0}$, interpolate the obtained data and sum on all values of
$x_{0}$ to have the back-projected data $\mathcal{G}^{\ddagger}f(x,z)$;

4 Perform the 2D FFT of $\mathcal{G}^{\ddagger}f(x,z)$ and weight by
$2\pi^{2}|\sigma|$;

5 Compute the inverse FFT of the result to recover $f$ ;

Algorithm 1 Reconstruction of object $f$

##  5 Simulations results

The original object used for simulations is Derenzo phantom, an object made of
multiple circles of different sizes. The circles in the object also allow the
study of the performance of the algorithm in front of different contrasts and
spatial resolution, as well as its ability to reconstruct features locally
tangent to lines of any slope. The unit length used here is the pixel. We
suppose thus that the distance between the source and the detector paths of
the modality is two pixels. The size of the object is $N\times N=256\times
256$ pixels in all simulations. Furthermore, given the linear geometry of this
modality, there is no loss of generality to consider the object centred
relatively to the $z$-axis.

###  5.1 Data acquisition

Data measurement is calculated according to (4.1). Figure 4 shows an example
of data obtained for Derenzo phantom.

(a)

(b)

Figure 4: (a) Original object: Derenzo phantom. (b) Corresponding acquired
data for $N_{SD}=2048$ and $N_{r}=1024$. A distance of one pixel is left
between the upper part of the image and the detector path ($\delta=1$, see
5.2.1).

###  5.2 Influence of some parameters on reconstruction quality

In the following paragraphs, we study the influence of the other general
parameters of the system such as the position of the object, the number of
required positions for sensors or the number of scanning circle arcs. In
addition to a visual comparison of the reconstruction quality, we propose here
to measure quantitatively the error rate between the original object $f_{0}$
and the reconstruction $f$ with the Normalized Mean Squared Error (NMSE)
$=||f-f_{0}||_{2}^{2}/N^{2}$ where $||.||_{2}$ refers to the $2$-norm.

The regularization parameter $\epsilon$, which has to be small, was
arbitrarily set to $0.01$.

####  5.2.1 Position of the object relative to the detector path

We analysed here the influence of the position of the object on reconstruction
quality. Firstly, the number of positions $N_{SD}$ for the pair source -
detector and the number of scanning circles per position $N_{r}$ was chosen to
largely satisfy the well-known condition [33] $N_{SD}\times N_{r}\geq N^{2}$
and are set arbitrarily to $N_{SD}=1024$ and $N_{r}=1024$. A convenient choice
for these parameters will be discussed later. We performed various
acquisition, modifying the gap $\delta$ between the detector path and the
upper part of the object (see Figure 3). Consequently, the object is in the
square of Cartesian coordinates
$\left(x\in\left[-\frac{N}{2}+1,\frac{N}{2}\right];z\in[-N-\delta,-\delta+1\right])$.
Figure 5 shows the reconstruction results for $\delta=1,26$ and $51$ pixels
which correspond respectively a position for the object in the respective
domains $[-N-1,0]$, $[-N-26,-25]$ and $[-N-51,-50]$ along the z-axis. The
difference between the three reconstructions is on the top of the object. If
the object is close to the line of movement of the detector, then this part is
less well reconstructed. Indeed, this distance between the object and the
detector path allows having arcs of circle tangent horizontally to this part
of the object. An offset of $\delta=51$ seems to be a good trade-off between
quality of reconstruction and the applicability of such a measure in practical
use. For the rest of the simulations, $\delta$ is set to $51$.

(a) Reconstruction for $\delta=1$
NMSE = 0.0112

(b) Reconstruction for $\delta=26$
NMSE = 0.0074

(c) Reconstruction for $\delta=51$
NMSE = 0.0061

Figure 5: Reconstruction results of the Derenzo phantom 4a for $\delta=1$ (a),
$\delta=26$ (b) and $\delta=51$ (c) pixel(s).

####  5.2.2 Number of necessary positions for the pair source-detector

We studied then the number of different positions required for a good quality
of reconstruction. The influence of two running parameters is analysed, first,
the farthest position $x_{0,max}$ from the object for the source-detector pair
(that is an array of length $[-x_{0,max},x_{0,max}]$ for the source and
detector paths) and the distance $\Delta_{x_{0}}$ between two adjacent
positions of the pair.

Figures 6a, 6b and 6c show the reconstruction results when the farthest
position from the object to the pair source detector is respectively $2N,3N$
and $4N$ with a common $\Delta_{x_{0}}$ set to $1$. Reconstruction from a
domain $[-x_{0,max},x_{0,max}]=[-2N,2N]$ appears to be blurred with strong
artefacts in the upper parts of the image. For $x_{0,max}=3N$ and $4N$,
reconstruction quality seems to be visually equivalent, even if the NMSE for
$x_{0,max}=4N$ is higher. This may be due to numerical approximations. For the
rest of the simulation, $x_{0,max}$ is set to $3N$.

The influence of the distance $\Delta_{x_{0}}$ between two adjacent positions
of the pair is now evaluated. Figures 6d, 6e and 6f show the result for
$\Delta_{x_{0}}=2,1$ and $0.5$ pixels. This represents a respective amount of
$0.5,1$ and $2$ detectors per unit length. In Figure 6d $(\Delta_{x_{0}}=2)$,
we can see streaks suggesting a lack of data for reconstructing the object. On
the contrary, the doubling of the number of detectors between Fig. 6e and Fig.
6f does not bring a better quality of reconstruction. Consequently, the use of
one detector per unit length seems to be a good trade-off, and the value will
remain constant in the rest of the paper.

(a) Reconstruction for $x_{0,max}=2N$. NMSE = 0.0084

(b) Reconstruction for $x_{0,max}=3N$. NMSE = 0.0049

(c) Reconstruction for $x_{0,max}=4N$. NMSE = 0.0061

(d) Reconstruction for $0.5$ detectors per unit length. NMSE = 0.0046

(e) Reconstruction for $1$ detectors per unit length. NMSE = 0.0049

(f) Reconstruction for $2$ detectors per unit length. NMSE = 0.0049

Figure 6: Evaluation of the number of source-detector positions on
reconstruction quality. First row: Reconstruction results of the Derenzo
phantom 4a for $x_{0,max}=2N$ 6a, $3N$ 6b and $4N$ 6c where
$\Delta_{x_{0}}=1$. Second row: Reconstruction results for $0.5$ 6d, $1$ 6e
and $0.5$ 6f detector per unit length and $x_{0,max}=3N$ remains constant.

####  5.2.3 Number of scanning circles per position of the pair source-
detector

The number of scanning circles necessary for reconstruction is now under
study. In the same way, two parameters are of interest, that is, the maximum
radius $r_{max}$ of the scanning double circle arcs to be taken and the
discretization step $\Delta_{r}$ that have to be chosen. We first evaluate the
consequences of the value of $r_{max}$ with three examples on Fig. 7a, 7b and
7c where $r_{max}$ is set respectively to $2N,3N$ and $4N$ and $\Delta_{r}=1$.
For $r_{max}=2N$, the reconstruction suggests a lack of data in front of the
obtained results for $r_{max}=3N$ and $r_{max}=4N$. Notice the higher NMSE for
$r_{max}=4N$, probably due to numerical approximations.

We were finally looking for the appropriate discretization step for $r$,
setting $\Delta r$ to $1,2$ and $4$. Reconstruction results are shown
respectively in Fig. 7d, 7e and 7f. Reconstructions with a large
discretization step exhibit blur.

(a) Reconstruction for $r_{max}=2N$. NMSE = 0.0058

(b) Reconstruction for $r_{max}=3N$. NMSE = 0.0040

(c) Reconstruction for $r_{max}=4N$. NMSE = 0.0049

(d) Reconstruction for $\Delta_{r}=1$. NMSE = 0.0040

(e) Reconstruction for $\Delta_{r}=2$. NMSE = 0.0043

(f) Reconstruction for $\Delta_{r}=4$. NMSE = 0.0043

Figure 7: Evaluation of the number of scanning circles on reconstruction
quality. First row: Reconstruction results of the Derenzo phantom 4a for
$r_{max}=2N$ 7a, $3N$ 7b and $4N$ 7c where $\Delta_{r}=1$. Second row:
Reconstruction results for $\Delta_{r}=1$ 7d, $2$ 7e and $4$ 7f detector per
unit length and $r_{max}=3N$ remains constant.

####  5.2.4 Discussions

The above simulations results show some interesting issues for the CST
modality and the reconstruction quality that can be expected with such a
system. First, this system is able to scan objects whose depth is largely
oversized relative to the source-detector distance. However, this seems to be
counterbalanced by the need for a consequent length for the source and
detector linear paths, since sufficient reconstruction results appear for
source and detector paths of length size six times greater than the object
size. The necessary number of scanning circle arcs is also very important
since, if we relate the values of $r_{max}$ with the scattering angles, the
reconstruction quality is largely improved when $r_{max}$ is high whereas the
angular distance between two values of $r_{max}$ is very small. Considering a
larger distance between the source and the detector paths for the system may
reduce partially the required amount of projection data.

Moreover, the geometry offers a sufficient reconstruction quality for every
tangent with arbitrary slopes. However, one can notice that vertical slopes
are slightly less well reconstructed than the other ones. This can be seen if
we pay carefully attention to the left and right sides of the reconstruction.
This may be problematic if the object to scan is essentially made of vertical
features. One way to avoid this is to perform additional scans by rotating the
object, if possible.

Some artefacts that look like shadows around the circles which compose the
object remains clearly visible. The issue of artefacts in CST has already been
addressed in different manners, for instance in [36], where microlocal
analysis was employed to alleviate artefacts in reconstructions from limited
data. Moreover, in [37], a penalized iterative algorithm was developed for a
mixed modality. Regarding our approach, it can be combined in a pipeline with
post-processing stages based on machine learning, as we did in [38] for
limited data issues in classical computed tomography. Some work is on the way.

##  6 Concluding remarks

In this article, we proposed a new reconstruction algorithm for a recently
proposed CST modality with translational geometry. The algorithm can be
numerically implemented efficiently and reconstructions exhibit good quality.
A quantitative study of the required data has also been carried out. This
study proved the ability for such a CST system to reconstruct larger objects
than the system itself. This advantage is moderated by the fact that it is
necessary to have an important amount of source-detector positions, which can
make the acquisition time longer. An interesting issue of this work concerns a
reconstruction algorithm for the three-dimensional extension of this system
with planar paths for the source and the detector.

##  7 Acknowledgements

We would like to thank Prof. T. T. Truong for stimulating discussions.

C. Tarpau research work is supported by grants from Région Île-de-France (in
Mathematics and Innovation) 2018-2021 and LabEx MME-DII (Modèles Mathématiques
et Économiques de la Dynamique, de l’Incertitude et des Interactions) (No.
ANR-11-LBX-0023-01).

J. Cebeiro research work is supported by a postdoctoral grant from the
University of San Martín. He is also partially supported by SOARD-AFOSR (grant
number FA9550-18-1-0523).

## References

  * [1] P. Lale, “The examination of internal tissues, using gamma-ray scatter with a possible extension to megavoltage radiography,” _Physics in Medicine & Biology_, vol. 4, no. 2, p. 159, 1959.
  * [2] R. Clarke and G. Van Dyk, “Compton-scattered gamma rays in diagnostic radiography,” in _Medical Radioisotope Scintigraphy. VI Proceedings of a Symposium on Medical Radioisotope Scintigraphy_ , 1969.
  * [3] F. Farmer and M. P. Collins, “A new approach to the determination of anatomical cross-sections of the body by compton scattering of gamma-rays,” _Physics in Medicine & Biology_, vol. 16, no. 4, p. 577, 1971.
  * [4] G. Redler, K. C. Jones, A. Templeton, D. Bernard, J. Turian, and J. C. Chu, “Compton scatter imaging: A promising modality for image guidance in lung stereotactic body radiation therapy,” _Medical physics_ , vol. 45, no. 3, pp. 1233–1240, 2018.
  * [5] K. C. Jones, G. Redler, A. Templeton, D. Bernard, J. V. Turian, and J. C. Chu, “Characterization of Compton-scatter imaging with an analytical simulation method,” _Physics in Medicine & Biology_, vol. 63, no. 2, p. 025016, 2018\.
  * [6] S. Gautam, F. Hopkins, R. Klinksiek, and I. Morgan, “Compton interaction tomography I. Feasibility studies for applications in earthquake engineering,” _IEEE Transactions on Nuclear Science_ , vol. 30, no. 2, pp. 1680–1684, 1983.
  * [7] P. G. Prado, M. K. Nguyen, L. Dumas, and S. X. Cohen, “Three-dimensional imaging of flat natural and cultural heritage objects by a Compton scattering modality,” _Journal of Electronic Imaging_ , vol. 26, no. 1, p. 011026, 2017.
  * [8] G. Harding and E. Harding, “Compton scatter imaging: A tool for historical exploration,” _Applied Radiation and Isotopes_ , vol. 68, no. 6, pp. 993–1005, 2010.
  * [9] E. M. Hussein, M. Desrosiers, and E. J. Waller, “On the use of radiation scattering for the detection of landmines,” _Radiation Physics and Chemistry_ , vol. 73, no. 1, pp. 7–19, 2005.
  * [10] P. E. Cruvinel and F. A. Balogun, “Compton scattering tomography for agricultural measurements,” _Engenharia Agricola_ , vol. 26, no. 1, pp. 151–160, 2006.
  * [11] T. T. Truong and M. K. Nguyen, “Recent developments on Compton scatter tomography: theory and numerical simulations,” in _Numerical Simulation-From Theory to Industry_. IntechOpen, 2012.
  * [12] M. K. Nguyen and T. T. Truong, _Imagerie par rayonnement gamma diffusé_. Hermès Science, 2006\.
  * [13] J. Radon, “Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten,” _Akad. Wiss._ , vol. 69, pp. 262–277, 1917.
  * [14] A. M. Cormack, “Representation of a function by its line integrals, with some radiological applications,” _Journal of Applied Physics_ , vol. 34, no. 9, pp. 2722–2727, 1963.
  * [15] S. J. Norton, “Compton scattering tomography,” _Journal of applied physics_ , vol. 76, no. 4, pp. 2007–2015, 1994.
  * [16] M. K. Nguyen and T. T. Truong, “Inversion of a new circular-arc Radon transform for Compton scattering tomography,” _Inverse Problems_ , vol. 26, no. 6, p. 065005, 2010.
  * [17] G. Rigaud, M. K. Nguyen, and A. K. Louis, “Novel numerical inversions of two circular-arc Radon transforms in Compton scattering tomography,” _Inverse Problems in Science and Engineering_ , vol. 20, no. 6, pp. 809–839, 2012.
  * [18] G. Rigaud, R. Régnier, M. K. Nguyen, and H. Zaidi, “Combined modalities of Compton scattering tomography,” _IEEE Transactions on Nuclear Science_ , vol. 60, no. 3, pp. 1570–1577, 2013.
  * [19] T. T. Truong and M. K. Nguyen, “Radon transforms on generalized Cormack’s curves and a new Compton scatter tomography modality,” _Inverse Problems_ , vol. 27, no. 12, p. 125001, 2011.
  * [20] C. Tarpau, J. Cebeiro, M. Morvidone, and M. K. Nguyen, “A new concept of Compton Scattering tomography and the development of the corresponding circular Radon transform,” _IEEE Transactions on Radiation and Plasma Medical Sciences_ , vol. (accepted for publication), 2019, [10.1109/TRPMS.2019.2943555].
  * [21] C. Tarpau and M. K. Nguyen, “Compton scattering imaging system with two scanning configurations,” _Journal of Electronic Imaging_ , vol. 29, no. 1, p. 013005, 2020.
  * [22] J. Cebeiro, M. K. Nguyen, M. Morvidone, and C. Tarpau, “An interior Compton Scatter Tomography,” in _25th IEEE Nuclear Science Symposium and Medical Imaging Conference 2018 (IEEE NSS/MIC’18)_ , Sydney, Australia, Nov. 2018\.
  * [23] G. Rigaud, “Compton Scattering Tomography: Feature Reconstruction and Rotation-Free Modality,” _SIAM Journal on Imaging Sciences_ , vol. 10, no. 4, pp. 2217–2249, 2017. [Online]. Available: https://doi.org/10.1137/17M1120105
  * [24] C. Tarpau, J. Cebeiro, M. K. Nguyen, G. Rollet, and M. A. Morvidone, “Analytic inversion of a Radon transform on double circular arcs with applications in Compton Scattering Tomography,” _IEEE Transactions on Computational Imaging_ , vol. 6, pp. 958–967, 2020.
  * [25] J. Webber and E. L. Miller, “Compton scattering tomography in translational geometries,” _Inverse Problems_ , vol. 36, no. 2, p. 025007, 2020.
  * [26] J. W. Webber and W. R. Lionheart, “Three dimensional Compton scattering tomography,” _Inverse Problems_ , vol. 34, no. 8, p. 084001, 2018.
  * [27] G. Rigaud and B. N. Hahn, “3D Compton scattering imaging and contour reconstruction for a class of Radon transforms,” _Inverse Problems_ , vol. 34, no. 7, p. 075004, 2018.
  * [28] J. Cebeiro, C. Tarpau, M. A. Morvidone, D. Rubio, and M. K. Nguyen, “On a three dimensional Compton scattering tomography system with fixed source,” _Inverse Problems_ , vol. 37, no. 5, p. 054001, 2021.
  * [29] J. W. Webber and S. Holman, “Microlocal analysis of a spindle transform,” _AIMS Inverse Problems and Imaging_ , vol. 13, no. 2, pp. 231–261, 2019. [Online]. Available: http://aimsciences.org//article/id/7ad5560c-e076-4384-9e9d-1dab4121da6d
  * [30] J. Cebeiro, M. K. Nguyen, M. Morvidone, and A. Noumowé, “New “improved” Compton scatter tomography modality for investigative imaging of one-sided large objects,” _Inverse Problems in Science and Engineering_ , vol. 25, no. 11, pp. 1676–1696, 2017.
  * [31] T. Truong and M. Nguyen, “Compton scatter tomography in annular domains,” _Inverse Problems_ , vol. 35, no. 5, p. 054005, 2019.
  * [32] I. S. Gradshteyn, I. M. Ryzhik, D. Zwillinger, and V. Moll, _Table of integrals, series, and products; 8th ed._ Amsterdam: Academic Press, Sep 2014. [Online]. Available: https://cds.cern.ch/record/1702455
  * [33] R. N. Bracewell, “Numerical transforms,” _Science_ , vol. 248, no. II May, pp. 697–704, 1990.
  * [34] H. Bateman, _Tables of Integral Transforms_. New York: McGraw-Hill Book Compagny, 1954, vol. 1.
  * [35] T. T. Truong, “Function reconstruction from reflection symmetric radon data,” _Symmetry_ , vol. 12, no. 6, 2020. [Online]. Available: https://www.mdpi.com/2073-8994/12/6/956
  * [36] J. W. Webber and E. T. Quinto, “Microlocal analysis of a compton tomography problem,” _SIAM Journal on Imaging Sciences_ , vol. 13, no. 2, pp. 746–774, 2020.
  * [37] J. W. Webber, E. T. Quinto, and E. L. Miller, “A joint reconstruction and lambda tomography regularization technique for energy-resolved x-ray imaging,” vol. 36, no. 7, p. 074002, jul 2020. [Online]. Available: https://doi.org/10.1088/1361-6420/ab8f82
  * [38] I. Ayad, C. Tarpau, M. K. Nguyen, and N. S. Vu, “Deep morphological network-based artifact suppression for limited-angle tomography,” in _Proceedings of the 25th International Conference on Image Processing, Computer Vision and Pattern Recognition (IPCV’21)_ , Las Vegas, United States, Jul. 2021.