From 4549204507a44b648785272cbae43d5221e8a28b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sun, 14 Jul 2013 17:39:48 +0200 Subject: [PATCH] iradon: Clean up filter code. --- skimage/transform/radon_transform.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index db43090b..fc827c8b 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -15,7 +15,7 @@ References: """ from __future__ import division import numpy as np -from scipy.fftpack import fftshift, fft, ifft +from scipy.fftpack import fft, ifft, fftfreq from scipy.interpolate import interp1d from ._warps_cy import _warp_fast from ._radon_transform import sart_projection_update @@ -203,26 +203,26 @@ def iradon(radon_image, theta=None, output_size=None, img = util.pad(radon_image, pad_width, mode='constant', constant_values=0) # Construct the Fourier filter - f = fftshift(abs(np.mgrid[-1:1:2 / projection_size_padded])).reshape(-1, 1) - w = 2 * np.pi * f - # Start from first element to avoid divide by zero + 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": - f[1:] = f[1:] * np.sin(w[1:] / 2) / (w[1:] / 2) + # Start from first element to avoid divide by zero + fourier_filter[1:] = fourier_filter[1:] * np.sin(omega[1:]) / omega[1:] elif filter == "cosine": - f[1:] = f[1:] * np.cos(w[1:] / 2) + fourier_filter *= np.cos(omega) elif filter == "hamming": - f[1:] = f[1:] * (0.54 + 0.46 * np.cos(w[1:])) + fourier_filter *= (0.54 + 0.46 * np.cos(omega / 2)) elif filter == "hann": - f[1:] = f[1:] * (1 + np.cos(w[1:])) / 2 + fourier_filter *= (1 + np.cos(omega / 2)) / 2 elif filter is None: f[1:] = 1 else: raise ValueError("Unknown filter: %s" % filter) - filter_ft = np.tile(f, (1, len(theta))) # Apply filter in Fourier domain - projection = fft(img, axis=0) * filter_ft + projection = fft(img, axis=0) * fourier_filter radon_filtered = np.real(ifft(projection, axis=0)) # Resize filtered image back to original size