scikit-image · kczimm · Nov 23, 2017 · Nov 23, 2017 · Nov 23, 2017 · Nov 23, 2017
diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py
@@ -1,7 +1,7 @@
 from .hough_transform import (hough_line, hough_line_peaks,
                               probabilistic_hough_line, hough_circle,
                               hough_circle_peaks, hough_ellipse)
-from .radon_transform import (radon, iradon, iradon_sart,
+from .radon_transform import (radon, iradon, iradon_filter, iradon_sart,
                               order_angles_golden_ratio)
 from .finite_radon_transform import frt2, ifrt2
 from .integral import integral_image, integrate
@@ -26,6 +26,7 @@
            'hough_line_peaks',
            'radon',
            'iradon',
+           'iradon_filter',
            'iradon_sart',
            'order_angles_golden_ratio',
            'frt2',

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
@@ -117,6 +117,67 @@ def _sinogram_circle_to_square(sinogram):
     return np.pad(sinogram, pad_width, mode='constant', constant_values=0)
 
 
+def iradon_filter(filter_name, size):
+    """
+    Create iradon filter.
+
+    Generate a reconstruction kernel to be used in iradon during
+    the filtered back-projection algorithm.
+
+    Parameters
+    ----------
+    filter_name : str
+        Name of the filter.
+    size : int
+        Length of the filter array.
+
+    Returns
+    -------
+    kernel : array_like, dtype=float
+        Reconstruction filter to be used in iradon.
+
+    References
+    ----------
+    .. [1] AC Kak, M Slaney, "Principles of Computerized Tomographic
+           Imaging", IEEE Press 1988.
+    """
+    # Construct the Fourier filter
+
+    # The ramp filter is implemented based on eq. 61 on  p.72 of Kak and Slaney
+    # https://engineering.purdue.edu/~malcolm/pct/CTI_Ch03.pdf
+    n1 = np.arange(0, size / 2 + 1, dtype=np.int)
+    n2 = np.arange(size / 2 - 1, 0, -1, dtype=np.int)
+    n = np.concatenate((n1, n2))
+    f = np.zeros(size)
+    f[0] = 0.25
+    f[1::2] = -1 / (np.pi * n[1::2])**2
+
+    omega = 2 * np.pi * fftfreq(size)
+    kernel = 2 * np.real(fft(f))         # ramp filter
+    if filter_name == "ramp":
+        pass
+    elif filter_name == "shepp-logan":
+        # Start from first element to avoid divide by zero
+        kernel[1:] *= np.sin(omega[1:] / 2) / (omega[1:] / 2)
+    elif filter_name == "cosine":
+        freq = (0.5 * np.arange(0, size)
+                / size)
+        cosine_filter = np.fft.fftshift(np.sin(2 * np.pi * np.abs(freq)))
+        kernel *= cosine_filter
+    elif filter_name == "hamming":
+        hamming_filter = np.fft.fftshift(np.hamming(size))
+        kernel *= hamming_filter
+    elif filter_name == "hann":
+        hanning_filter = np.fft.fftshift(np.hanning(size))
+        kernel *= hanning_filter
+    elif filter_name is None:
+        kernel[:] = 1
+    else:
+        raise ValueError("Unknown filter: %s" % filter_name)
+
+    return kernel[:, np.newaxis]
+
+
 def iradon(radon_image, theta=None, output_size=None,
            filter="ramp", interpolation="linear", circle=True):
     """
@@ -203,25 +264,8 @@ def iradon(radon_image, theta=None, output_size=None,
     pad_width = ((0, projection_size_padded - radon_image.shape[0]), (0, 0))
     img = np.pad(radon_image, pad_width, mode='constant', constant_values=0)
 
-    # Construct the Fourier filter
-    f = fftfreq(projection_size_padded).reshape(-1, 1)   # digital frequency
-    omega = 2 * np.pi * f                                # angular frequency
-    fourier_filter = 2 * np.abs(f)                       # ramp filter
-    if filter == "ramp":
-        pass
-    elif filter == "shepp-logan":
-        # Start from first element to avoid divide by zero
-        fourier_filter[1:] = fourier_filter[1:] * np.sin(omega[1:]) / omega[1:]
-    elif filter == "cosine":
-        fourier_filter *= np.cos(omega)
-    elif filter == "hamming":
-        fourier_filter *= (0.54 + 0.46 * np.cos(omega / 2))
-    elif filter == "hann":
-        fourier_filter *= (1 + np.cos(omega / 2)) / 2
-    elif filter is None:
-        fourier_filter[:] = 1
-    else:
-        raise ValueError("Unknown filter: %s" % filter)
+    fourier_filter = iradon_filter(filter, projection_size_padded)
+
     # Apply filter in Fourier domain
     projection = fft(img, axis=0) * fourier_filter
     radon_filtered = np.real(ifft(projection, axis=0))

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
@@ -4,7 +4,7 @@
 import numpy as np
 from skimage import data_dir
 from skimage.io import imread
-from skimage.transform import radon, iradon, iradon_sart, rescale
+from skimage.transform import radon, iradon, iradon_sart, rescale, iradon_filter
 
 from skimage._shared import testing
 from skimage._shared.testing import test_parallel
@@ -45,6 +45,27 @@ def _rescale_intensity(x):
     return x
 
 
+def test_iradon_rampfilter_bias_circular_phantom():
+    """
+    test that a uniform circular phantom has a small reconstruction bias using
+    the ramp filter
+    """
+    pixels = 128
+    xy = np.arange(-pixels / 2, pixels / 2) + 0.5
+    x, y = np.meshgrid(xy, xy)
+    image = x**2 + y**2 <= (pixels/4)**2
+
+    theta = np.linspace(0., 180., max(image.shape), endpoint=False)
+    sinogram = radon(image, theta=theta)
+
+    reconstruction_fbp = iradon(sinogram, theta=theta)
+    error = reconstruction_fbp - image
+
+    tol = 5e-5
+    roi_err = np.abs(np.mean(error))
+    assert( roi_err < tol )
+
+
 def check_radon_center(shape, circle):
     # Create a test image with only a single non-zero pixel at the origin
     image = np.zeros(shape, dtype=np.float)