# Exact sinogram: an analytical approach
to the Radon transform of phantoms

Monica Dessole 1 1 Department of Mathematics “Tullio Levi Civita”, University
of Padova, Italy.
2 PNC - Padova Neuroscience Center, Italy.  Marta Gatto 1 1 Department of
Mathematics “Tullio Levi Civita”, University of Padova, Italy.
2 PNC - Padova Neuroscience Center, Italy.  Davide Poggiali 2 1 Department of
Mathematics “Tullio Levi Civita”, University of Padova, Italy.
2 PNC - Padova Neuroscience Center, Italy.  Francesca Tedeschi 1 1 Department
of Mathematics “Tullio Levi Civita”, University of Padova, Italy.
2 PNC - Padova Neuroscience Center, Italy.

###### Abstract

Phantoms can serve as a gold standard for the validation of MRI numerical
methods. In some special cases, it is possible to compute analytically the
Radon transform, or sinogram, of a phantom. In this work, we present
analytical formulae to compute the exact sinograms of three classes of
phantoms. We compare the use of the discrete Radon transform, that yields an
approximate sinogram, and the correspondent analytical sinogram for image
reconstruction.

##  1 Introduction and motivations

The imaging technique of exposing an object with rays at different angles is
called _tomography_ , and it produces a projection image of the inaccessible
regions of the body. The regions are determined by their attenuation function,
that quantifies how the intensity of radiations is attenuated, proportionally
to the density, due to mass absorption. The inverse problem of recovering a
function from its integral values along cross sections was first studied by
Johann Radon in 1917 [7], implemented by Cormack [1] to Computerized
Tomography (CT) technologies in the 1960s. The graphic representation of the
Radon transform is commonly referred to as _sinogram_.
Numerical methods for image reconstruction are tested on the so-called
_phantoms_ , fictitious images that reproduce a body section, as the Shepp-
Logan phantom [8]. Using a phantom gives the advantage to test a
reconstruction algorithm on a zero-noise data so the error we get is only due
to numerical inaccuracies in the algorithm itself [3, 6, 4]. In some special
cases, it is possible to compute analytically the Radon transform of a
phantom.
The scope of this work is to compare the use of the discrete Radon transform,
that yields an approximate sinogram, and the correspondent analytical sinogram
for image reconstruction.
We created exact-sinogram [2], an open-source python library of Analytical
Radon Transforms for different classes of phantoms, obtained as a combination
of ellipses, squares or rectangles, in order to successively reconstruct the
original image and to analyze the error.

##  2 Mathematical tools

Let $f=f(x,y)$ an integrable continuous function with compact support,
standing for the attenuation function, with
$(x,y)\in\Omega\subset\mathbb{R}^{2}$. We define the Radon transform of $f$,
as the following functional

$Rf(t,\theta):=\int_{l_{t,\theta}}f\,ds=\int_{s=-\infty}^{\infty}f(t\cos{\theta}-s\sin{\theta},t\sin{\theta}+s\cos{\theta})\,ds.\qquad\qquad\forall\,(t,\theta)\in\mathbb{R}\times[0,2\pi),$

where $l_{t,\theta}$ is the line that passes through the point
$(t\cos(\theta),t\sin(\theta))$ and orthogonal to to the unit vector
$\mathbf{n}=(\cos(\theta),\sin(\theta))$. The classical way of representing
the Radon transform is by treating $(t,\theta)$ as rectangular coordinates.
The values of the Radon transform are depicted in a continuum grey scale from
$0$ to $1$, standing for black and white respectively. This picture is the
sinogram and it represents data generated by the X-ray emission/detection
machine for the given slice of the sample.

We aim to reconstruct an approximation of the function $f$ from a set of
values of the Radon transform calculated on a discrete grid $(t,\theta)$,
using mathematical phantoms. A phantom is an image in grey scale with colors
representing mass absorption, made of a collection of geometric figures
(possibly convex), each of which has a different constant attenuation.

In this paper we exploit the analytical Radon transform of two types of
figures: the ones that can be described by a single mathematical equation,
e.g. ellipses, and the polygons, e.g. squares and rectangles.

We give here a brief description of two examples. More details can be found in
[3].

######  Phantom (Phantom of ellipses).

Every ellipse in the $xy$-plane is uniquely determined by six real parameters:
its semi-axes $a$ and $b$, the coordinates ($x_{0},y_{0}$) of its center
point, the angle $\phi$ of rotation with respect to the $x$-axis, and the
value $\delta$ of the attenuation function in the region of the plane inside
the ellipse. The coordinates in the roto-translated frame are
$\hat{x}=(x-x_{0})\cos{\phi}+(y-y_{0})\sin{\phi}$ and
$\hat{y}=(y-y_{0})\cos{\phi}-(x-x_{0})\sin{\phi}$.
Now, let $\epsilon$ be a generic ellipse and for real numbers $t$ and
$\theta$, let pose $\hat{\theta}=\theta-\phi$ and
$\hat{t}=t-x_{0}\cos{\theta}-y_{0}\sin{\theta}$. We assume the attenuation
function has the following form

$f_{\epsilon}(x,y)\coloneqq\begin{cases}\delta&\text{if
}\frac{\hat{x}^{2}}{a^{2}}+\frac{\hat{y}^{2}}{b^{2}}\leq 1,\\\
0&\text{otherwise}.\end{cases}$

The Radon transform of $f_{\epsilon}$ is obtained computing analytically the
intersections between the line $l_{t,\theta}$ and the ellipse, making the
difference between them, i.e. their distance, and is expressed as

$\begin{cases}Rf_{\epsilon}(t,\theta)=\delta\frac{2ab\sqrt{b^{2}\sin^{2}{\hat{\theta}}+a^{2}\cos^{2}{\hat{\theta}}-\hat{t}^{2}}}{b^{2}\sin^{2}{\hat{\theta}}+a^{2}\cos^{2}{\hat{\theta}}}&\text{if
}\hat{t}^{2}\leq b^{2}\sin^{2}{\hat{\theta}}+a^{2}\cos^{2}{\hat{\theta}},\\\
Rf_{\epsilon}(t,\theta)=0&\text{otherwise}.\end{cases}$ (1)

######  Phantom (Phantom of squares).

A square in general position in the plane is, instead, defined by five
parameters: coordinates $(x_{0},y_{0})$ of the center, side length $w$, angle
of counterclockwise rotation $\phi$ from the horizontal and the attenuation
value $\delta$. To generate the sinogram, we imagine the region inside the
square as the intersection of two infinite bands: one vertical and one
horizontal. Then, we compute the points at which the line crosses these bands
to find, among them, the “entry” and “exit” of the line within the square. The
distance between the exit and the entry points will be the Radon value
associated with those $(t,\theta)$, analogously to the case of the ellipse.

When two planar figures overlap in a phantom, the attenuation coefficient on
the overlap will be the sum of the individual attenuation coefficients. The
linearity property of the Radon transform will allow us to compute the
analytical sinogram of any combination of such figures.

For each phantom, we are going to calculate two reconstructed images, one from
the approximated sinogram, computed with the discrete Radon transform, and
another from the exact sinogram. To do this, we apply the inverse Radon
transform provided by the iradon function of the Python library Scikit-image
[10]. We would expect a smaller error with the analytical approach.

Figure 1: Phantom of rectangles (left), its approximated (center) and exact
(right) sinograms.

##  3 The library

The python library exact-sinogram is available in Github
(https://github.com/francescat93/Exact_sinogram) or can be installed via pip.
It is written on top of numpy [9], scipy [11], scikit-image [10], and
matplotlib for displaying the images.

A phantom is an instance of the class `Phantom`, initialized with the phantom
type (`‘ellipses’`, `‘squares’`, `‘rectangles’`) and a matrix containing the
parameters of the planar figures to be displayed. We also provide a gallery of
default phantoms of each type, accessible by calling the function
`my_phantomgallery`. The class `Phantom` has two methods: `get_phantom` that
returns the image of the phantom, and `get_sinogram` generates the analytical
Radon transform, using a different procedure according to the phantom type.
The implementation of the transforms is provided by the Scikit-image library
[10].

The following lines of code provide an usage example of the exact-sinogram
package.

⬇

1import numpy as np

2from skimage.transform import radon, iradon

3

4import exact_sinogram as es

5

6circle = True

7phantom_type = ’ellipses’

8# phantom_type : ’ellipses’ (or ’shepp_logan’) , ’modified_shepp_logan’,
’squares’, ’rectangles’

9n_points = 300; # number of pixels

10

11# Creation of an instance of the Phantom Class

12Phm = es.Phantom(phantom_type = phantom_type, circle = circle)

13

14# Creation of the matrix-image of the phantom, w/ input the number of pixel:

15P = Phm.get_phantom(N = n_points)

16

17# Array of projection angles

18theta_vec_deg = np.linspace(0, 359, 360)

19

20# Exact sinogram and inversion (exact_sinogram package)

21## Calculate the excact Sinogram

22analytical_sinogram = Phm.get_sinogram(N = n_points, theta_vec =
np.deg2rad(theta_vec_deg))

23

24## invert the Radon transform on the exact sinogram

25P_an = iradon(analytical_sinogram, theta = theta_vec_deg, circle=circle)

26

27# Approximated sinogram and inversion (scikit-image)

28## Calculate the approximated sinogram

29approx_sinogram = radon(analytical_sinogram, theta = theta_vec_deg,
circle=circle)

30## invert the Radon transform on the approximated sinogram

31P_approx = iradon(approx_sinogram, theta = theta_vec_deg, circle=circle)

##  4 Results and conclusions

For each phantom type, we compute the error between the synthetic data and the
reconstructed data obtained by the following procedures:

  1. 1.

by computing the analytic sinogram and then applying the inverse Radon
transform;

  2. 2.

in the classic way, by using the Radon transform and its inverse.

We are going to evaluate the quality of the reconstructed image by measuring
the classical $2$-norm and we collect the results in Table 1.

| Ellipses (Shepp-Logan) | Ellipses (Modified Shepp-Logan) | Squares | Rectangles
---|---|---|---|---
$\frac{\|P_{an}-P\|_{2}}{\|P\|_{2}}$ | 0.03073 | 0.02590 | 0.02703 | 0.01838
$\frac{\|P_{inv}-P\|_{2}}{\|P\|_{2}}$ | 0.03604 | 0.03191 | 0.01997 | 0.02588
Table 1: Relative errors of the image reconstructed from the analytical
sinogram $P_{an}$ and the image reconstructed form the numerical sinogram
$P_{inv}$ with respect to the original image $P$, taking into account Gibbs
phenomenon.

The Gibbs phenomenon [5] is the concentration of the errors along the
discontinuities of the figure. With a mask, we eliminated the influence of
boundaries and appreciated our algorithm that gives a smaller relative error
(see Table 1) with respect to the reconstruction through the iradon function.

## References

  * Cormack [1963] A. M. Cormack.  Representation of a function by its line integrals, with some radiological applications.  _Journal of Applied Physics (U.S.)_ , 34, 9 1963.  doi: 10.1063/1.1729798.
  * Dessole et al. [2020] M. Dessole, M. Gatto, D. Poggiali, and F. Tedeschi.  Exact Sinogram.  https://github.com/francescat93/Exact_sinogram, 2020.
  * Feeman [2016] T. G. Feeman.  _The Mathematics of Medical Imaging_.  UNITEXT 92. Springer, 1 edition, 2016.
  * Hansen [2010] P. C. Hansen.  _Discrete Inverse Problems: Insight and Algorithms_.  Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2010.  ISBN 0898716969, 9780898716962.
  * Jerri [1998] A. J. Jerri.  _The Gibbs Phenomenon in Fourier Analysis, Splines and Wavelet Approximations_.  Mathematics and Its Applications 446. Springer US, 1 edition, 1998.  ISBN 978-1-4419-4800-7,978-1-4757-2847-7.
  * Kaipio and Somersalo [2005] J. Kaipio and E. Somersalo.  _Statistical and Computational Inverse Problems_.  Springer, Dordrecht, 2005.
  * Radon [1986] J. Radon.  On the determination of functions from their integral values along certain manifolds.  _IEEE Transactions on Medical Imaging_ , 5(4):170–176, 1986.
  * Shepp and Logan [1974] L. A. Shepp and B. F. Logan.  The Fourier reconstruction of a head section.  _IEEE Transactions on Nuclear Science_ , 21(3):21–43, 1974.
  * Van Der Walt et al. [2011] S. Van Der Walt, S. C. Colbert, and G. Varoquaux.  The numpy array: a structure for efficient numerical computation.  _Computing in Science & Engineering_, 13(2):22–30, 2011.  doi: 10.1109/MCSE.2011.37.
  * van der Walt et al. [2014] S. van der Walt, J. L. Schönberger, J. Nunez-Iglesias, F. Boulogne, J. D. Warner, N. Yager, E. Gouillart, and T. a. Yu.  scikit-image: image processing in python.  _PeerJ_ , 2:e453, June 2014.  ISSN 2167-8359.  doi: 10.7717/peerj.453.  URL https://doi.org/10.7717/peerj.453.
  * Virtanen et al. [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. Jarrod Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. Carey, İ. Polat, Y. Feng, E. W. Moore, J. Vand erPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and S. . . Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17:261–272, 2020.  doi: https://doi.org/10.1038/s41592-019-0686-2.