Fast(ish) approximation to a radon transform.

radon_transform.py

import numpy
import scipy.ndimage.interpolation

def radon_transform(ims):
    """Radon transform a batch of images.
    
    Parameters
    ----------
    ims : numpy.ndarray
      (N x W x H) NumPy array of N images.
    
    Returns
    -------
    (N x W x W) NumPy array of radon-transformed images.
    """
    projections = numpy.linspace(0, 360, ims.shape[1])
    rotations = numpy.transpose(numpy.stack([
        scipy.ndimage.interpolation.rotate(ims, projection, axes=(1, 2), reshape=False)
        for projection in projections]), (1, 0, 2, 3))
    transformed = rotations.sum(axis=3)
    return transformed