1
$\begingroup$

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.)

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
10
  • $\begingroup$ So what abut the answer you gave yourself on the thread you posted several years ago? Is it a good solution for find the Radon? $\endgroup$ – Ian Dec 3 '14 at 17:07
  • $\begingroup$ For that thread and the answer, the Matlab code you gave, what is (ceil(-M/2) : floor(M/2)-1) and (ceil(-M/2) : floor(M/2))? Why in parenthesis? There is no Matlab expressions like that... What is the imggrid function exactly? I still don't understand. Thanks! $\endgroup$ – Ian Dec 3 '14 at 17:11
  • $\begingroup$ @Ian It worked from what I remember. I'm not sure what the confusion about the matlab syntax is; it's perfectly valid. However, I added the contents of imggrid to my answer in that post. Hopefully that will clear things up. $\endgroup$ – AnonSubmitter85 Dec 3 '14 at 17:40
  • $\begingroup$ Does radon matrix just have value 0 and 1? Like the output of your code. Why? $\endgroup$ – Ian Dec 3 '14 at 17:45
  • $\begingroup$ @Ian Because all it is doing is summing a certain set of pixels together. Each output pixel in the radon transform is the sum of the pixels that lie along a certain line. To calculate this sum, we weight the pixels along this line by $1$ and those that don't by $0$, then we sum the weighted pixels together. $\endgroup$ – AnonSubmitter85 Dec 3 '14 at 17:48
1
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

Currently I am solving the same problem as you guys. I found a way to find the matrix.

First, we may view image(MxN) as a vector by certain reshape. Second, set the ONLY one entry to 1, and others to 0. Third, input this image to the function 'radon' then you way get the transformed image(LxP)

Loop over these three steps for every entry of image, which may causes MxN times. cascade the result image(LxP) MxN times, then you will get the answer.

I.e. for 2x3-->3x4
[1 0 0 0 0 0]-->X_1 (reshape to a row vector)
[0 1 0 0 0 0]-->X_2 (reshape to a row vector)
[0 0 1 0 0 0]-->X_3 (reshape to a row vector)
[0 0 0 1 0 0]-->X_4 (reshape to a row vector)
[0 0 0 0 1 0]-->X_5 (reshape to a row vector)
[0 0 0 0 0 1]-->X_6 (reshape to a row vector)

Finally, stack result row vectors
[X_1, X_2, X_3, X_4, X_5, X_6].T
$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.