PDF Print

Descrete Base for the reconstruction

Let us take a linear operator that describes the image formation. Let our starting space be spanned by the n-dimensional real numbers and this operator should transform functions with variables in this space into the acquisition data sinogram space. For example, and for showing the underlying idea we could think that the Radon transform is such a linear operator. Let us define a complete orthogonal base with elements $ \varphi _{i}. Now for a sufficiently well behaving f(x) function:
$  f=\sum_{i}c_{i}\varphi _{i} (1)
If our operator is linear:
$ \mathfrak{R}f=\sum_{i}c_{i}\mathfrak{R}\varphi _{i}(2)
In physics it is common to use series expansion for solving differential equations. Physical quantities are often continuous, therefore e.g. the Fourier expansion or polynomial expansion (Legendre, Chebyshev etc.) or their combination (e.g. spherical harmonics) would be applied. There are successful attempts at approximate inverse Radon transforms using spherical harmonics but it has not gained much popularity up till now.

If we have an orthonormal base, then by definition:
$ \int \varphi _{i}\varphi _{j}=\delta _{ij} (3)
multiplyin both sides of Eq. (2)
$\varphi _{i}The coefficient can be calculated as follows:
$ c_{i}=\int \varphi _{i} f (4)

Eq. (3) and Eq. (4) are not only valid for orthogonal function bases but we could also take a discretisation base $\Pi\left ( \mathbf{x} \right )_{i}, containing characteristic functions that are in the "i"th domain
$1/\sqrt{\left | \Delta x_{i}  \right |}, and otherwise 0; here $\left | \Delta x_{i}  \right | the measure of the i. domain. Eq. (3) and Eq. (4) holds if we restrict the function base to be disjoint. For easier calculations a structured mesh is used, like a simple cubic grid. If f continuous, then Eq. (1) is only valid approxiomately. Now, for N number of domains:

$ f=\sum_{i}^{N}\Pi_{i}\left ( \mathbf{x} \right )c_{i}

Note, that with this base the expansion of f means the integral average taken on the ith domain.

This base functions are called volumetric pixels or in short: voxels. If talked about 3D we could use 3 indices according to the 3 spatial coordinates. It is possible to count all of them using a single index, and this is the notation we are going to use, where the i index runs through every domain in every dimension.

We can use the same discretisation in the sinogram space spanned by $  t,\boldsymbol{\omega } . Let us consider this notation in a general sense, i.e. let it also be meant for the variables of the ray transform. On this space we can use the notation for this expansion like this:
$\Pi_{j}\left ( t,\boldsymbol{\omega } \right ).
In practice the expansion of the sinogram space means actually the expansion of the ray-transform, e.g. for x-ray systems the line connecting the source and the detector pixel or in positron emission tomography the line connecting the crystals of a detector pair, called the Line of Response (LoR).

We can now write:
$ \mathfrak{R}f=\sum_{j}^{M}\Pi_{j}\left ( t,\boldsymbol{\omega } \right )y_{j}

The coefficients yj can be expressed as
$ y_{j}=\int \Pi_{j}\left ( t,\boldsymbol{\omega } \right )\mathfrak{R}f=\int \Pi_{j}\left ( t,\boldsymbol{\omega } \right )\mathfrak{R}\sum_{i}^{N}c_{i}\Pi_{i}\left (\mathbf{x}\right )=\sum_{i}^{N}c_{i}\underbrace{\int \Pi_{j}\left ( t,\boldsymbol{\omega } \right )\mathfrak{R}\Pi_{i}\left (\mathbf{x}\right )}_{A_{ij}}

It is also common to use instead of the ci coefficients the letter xi,indicating that these coefficients are the unknowns of the equation. With matrix notations our system becomes:

By the derivation we did not set the condition that the operator in question is the Radon transform, it is enough that it is a linear operator and projects into the sinogram space. Important aspect of the A mátrix that it is not square therefore cannot be inverted directly. In practice it is a sparse matrix, since a voxel only gives contribution to some LoRs that pass through that voxel. Its typical size e.g. for a pET scanner can be around 106x105.

In the next section we will discuss the matrix inversion solution to the above equation.

Site Language: English

Log in as…