What is the matrix representation of Radon transform?

Question

Just as the title, my question is what is the matrix representation of Radon transform (Radon projection matrix)? I want to have an exact matrix for the Radon transformation.

(I want to implement some electron tomography algorithms by myself so I need to use the matrix representation of the Radon transformation.)

(An introduction of the Radon transformation: click here.)

Answer 1

A couple of years ago I was looking for the same answer. You can see that thread here. However, I don't think my solution is the only way to do it nor necessarily the best way.

We know it's a linear transform, so if you can take the transform of the right identity matrix, you'll get the matrix representation. As an easy example using matlab notation, if we wanted to get the matrix representation of a column-wise DFT matrix using the fft() function, we could do it with A = fft(eye(N)). Now, it's a little more complicated with the radon transform to do it efficiently, but the same idea applies.

If the image is $\mathbf{X} \in \mathbb{C}^{M \times N}$ and we wanted to find the DRT matrix $\mathbf{R} : \mathbb{C}^{M \times N} \rightarrow \mathbb{C}^{L\times P}$, then we know that $\mathbf{R} \in \mathbb{C}^{LP \times MN}$ and that the DRT is given by $\mathbf{Y} = \mathbf{R} \cdot \operatorname{vec}(\mathbf{X}) \in \mathbb{C}^{LP \times 1}$, where $\operatorname{vec} : \mathbb{C}^{A\times B} \rightarrow \mathbb{C}^{AB \times 1}$. The typical 2-D display of the DRT would then be had by simply reshaping the $LP \times 1$ vector into a $L \times P$ array.

We don't know $\mathbf{R}$, but we have access to a function that will efficiently compute the DRT, call it radon(), such that Y = radon(X) is equivalent to $\mathbf{Y} = \mathbf{R} \cdot \operatorname{vec}(\mathbf{X})$. Since $\mathbf{R} \cdot \mathbf{I} \cdot \operatorname{vec}(\mathbf{X}) = \mathbf{R} \cdot \operatorname{vec}(\mathbf{X})$, where $\mathbf{I}$ is the identity matrix, if we can used radon() in such a way that A = radon(I), then we have that A * vec(X) = radon(X).

The problem is that I don't see how to use radon to efficiently compute the above. Most radon() implementations take a 2-D input and output a 2-D array. It's simple enough to reshape the input/output into vectors, but the problem is that you would have to do that $MN$ times. Even for modest image sizes, the above will become computationaly prohibitive rather fast.

I hope the above helps. If you come up with an efficient solution, please post it here or at the stackoverlow thread linked above.

Answer 2

Radon transform consists in a rotation plus a projection in X axis.

This code will throw True for any squared matrix where there is an inscribed circle containing the target values.

from skimage.transform import radon, rotate
import numpy as np


def get_image_center(img):

    center = (img.shape[0] * 0.5, img.shape[1] * 0.5)

    return center

def rotate_image(img, deg):

    rotated_img = rotate(img, 
                         deg, 
                         center=get_image_center(img), 
                         preserve_range=True)

    return rotated_img

theta_test = np.linspace(0, 360, 100)

sinogram = radon(image, theta=theta_test, circle=True)

def perform_radon(img, angle_vector):

    radon_result_list = []
    x_projection_matrix = np.ones(img.shape[0])

    for angle in angle_vector:
        rotated_img = rotate_image(img,  - angle)
        radon_transform = np.dot(x_projection_matrix, rotated_img)
        radon_result_list.append(radon_transform)

    arr = np.vstack(radon_result_list).T

    return arr

manual_radon = perform_radon(img, theta_test)
manual_radon.shape

np.allclose(manual_radon, sinogram)

The last thing needed to get a linear transformation that represents Radon transformation is to get a linear operator that represents the rotation.