The Radon transform is defined as
$$f(\theta,s) = \int_{-\infty}^\infty u\left(x(s,t), y(s,t)\right) \mathrm{d}t,$$with
$$x(s,t) = t \sin \theta + s \cos \theta, \quad y(s,t) = t\cos\theta - s \sin \theta.$$We can now easily compute the matrix elements for a given set of angles $\{\theta_k\}_{k=1}^{n_\theta}$ and define
$$A = \left( \begin{matrix} A(\theta_1) \\ A(\theta_2) \\ \vdots \\ A(\theta_k)\end{matrix}\right).$$The vectorised discretised sinogram and image are then related via
$$\mathbf{f} = A\mathbf{u}.$$import numpy as np
x = np.array([-1, -.5, 0, .5, 1])
y = np.array([-1, -.5, 0, .5, 1])
h = 0.5
s = np.array([-1.25, -.75, -.25, .25, .75, 1.25])
theta = 1
si = -.5
tx = (x - si*np.cos(theta)) / np.sin(theta)
ty = -(y - si*np.sin(theta)) / np.cos(theta)
t = np.sort(np.concatenate((tx,ty)))
xi = si*np.cos(theta) + t*np.sin(theta)
yi = si*np.sin(theta) - t*np.cos(theta)
xc = h*((si*np.cos(theta) + 0.5*(t[0:-1]+t[1:])*np.sin(theta))//h) + h/2
yc = h*((si*np.sin(theta) - 0.5*(t[0:-1]+t[1:])*np.cos(theta))//h) + h/2
coo_array
Coordinate formatcsc_array
Compressed Sparse Column formatcsr_array
Compressed Sparse Row formatConsider a finite-difference matrix
$$A = \left(\begin{matrix}-2 & 1& 0 & \ldots & \\ 1 & -2 & 1 & 0 \\ & \ddots & \ddots & \ddots & \\ 0 & \ldots &0 & 1 & -2 \end{matrix}\right).$$The function would look like
def FD(u, n):
v = np.zeros(n)
v[0] = -2*u[0] + u[1]
v[1:-1] = u[:-2] - 2*u[1:-1] + u[2:]
v[-1] = u[-2] - 2*u[-1]
Consider summing the elements of a vector
$$A = (1, 1, 1, \ldots, 1).$$def sum(u,n):
return np.sum(u)
def sum_transpose(v,n):
u = v * np.ones(n)
Write a function RadonMatrix(n, s, theta)
that discretizes the Radon transform of a given n
$\times$ n
image, defined on $[-1,1]^2$, for given arrays s
and theta
and returns a matrix of size len(s)*len(theta)
by n*n
.
You can start with a loop that traces each ray individually and computes the weights in the matrix according to the procedure we outlined in class. If you are familiar with numpy
you can attempt to make it more efficient.
Write a function RadonSparseMatrix
which returns a sparse matrix representation of RadonMatrix
.
RadonMatrix
Write a function Radon(u, n, s, theta)
which performs the Radon transform of the given image u
and the corresponding transpose RadonTranspose(f, n, s, theta)
.