From 237fb989b0a76e6c67aa03d97952ef237527589c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Fri, 5 Jul 2013 14:05:06 +0200 Subject: [PATCH 1/8] radon: Use util.pad for array padding. --- skimage/transform/radon_transform.py | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index 944b4a6b..dcb13cae 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -18,6 +18,8 @@ import numpy as np from scipy.fftpack import fftshift, fft, ifft from ._warps_cy import _warp_fast from ._radon_transform import sart_projection_update +from .. import util + __all__ = ["radon", "iradon", "iradon_sart"] @@ -77,20 +79,14 @@ def radon(image, theta=None, circle=False): dh = padded_image.shape[0] // 2 dw = padded_image.shape[1] // 2 else: - height, width = image.shape diagonal = np.sqrt(2) * max(image.shape) - heightpad = int(np.ceil(diagonal - height)) - widthpad = int(np.ceil(diagonal - width)) - padded_image = np.zeros((int(height + heightpad), - int(width + widthpad))) - y0 = heightpad // 2 - y1 = y0 + height - x0 = widthpad // 2 - x1 = x0 + width - padded_image[y0:y1, x0:x1] = image + pad = [int(np.ceil(diagonal - s)) for s in image.shape] + pad_width = [(p // 2, p - p // 2) for p in pad] + padded_image = util.pad(image, pad_width, mode='constant', + constant_values=0) out = np.zeros((max(padded_image.shape), len(theta))) - dh = y0 + height // 2 - dw = x0 + width // 2 + dh = pad[0] // 2 + image.shape[0] // 2 + dw = pad[1] // 2 + image.shape[1] // 2 shift0 = np.array([[1, 0, -dw], [0, 1, -dh], From ac460a9777303df3d4df17a06c1ec70b53122f68 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Fri, 5 Jul 2013 14:13:41 +0200 Subject: [PATCH 2/8] iradon: Use util.pad for array padding. This depends on util.pad accepting 0s in pad_width, which has a proposed solution in gh-634. --- skimage/transform/radon_transform.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index dcb13cae..93979977 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -109,11 +109,10 @@ def radon(image, theta=None, circle=False): def _sinogram_circle_to_square(sinogram): - size = int(np.ceil(np.sqrt(2) * sinogram.shape[0])) - sinogram_padded = np.zeros((size, sinogram.shape[1])) - pad = (size - sinogram.shape[0]) // 2 - sinogram_padded[pad:pad + sinogram.shape[0], :] = sinogram - return sinogram_padded + diagonal = int(np.ceil(np.sqrt(2) * sinogram.shape[0])) + pad = diagonal - sinogram.shape[0] + pad_width = ((pad // 2, pad - pad // 2), (0, 0)) + return util.pad(sinogram, pad_width, mode='constant', constant_values=0) def iradon(radon_image, theta=None, output_size=None, From 5955e4e612b7efdd3eafc0f105855d13959fff5d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Fri, 5 Jul 2013 15:27:39 +0200 Subject: [PATCH 3/8] iradon: Reduce code duplication. --- skimage/transform/radon_transform.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index 93979977..71548324 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -232,8 +232,6 @@ def iradon(radon_image, theta=None, output_size=None, k = np.round(mid_index + ypr * np.cos(th[i]) - xpr * np.sin(th[i])) backprojected = radon_filtered[ ((((k > 0) & (k < n)) * k) - 1).astype(np.int), i] - if circle: - backprojected[~reconstruction_circle] = 0. reconstructed += backprojected elif interpolation == "linear": for i in range(len(theta)): @@ -244,11 +242,11 @@ def iradon(radon_image, theta=None, output_size=None, b1 = ((((b > 0) & (b < n)) * b) - 1).astype(np.int) backprojected = (t - a) * radon_filtered[b0, i] + \ (a - t + 1) * radon_filtered[b1, i] - if circle: - backprojected[~reconstruction_circle] = 0. reconstructed += backprojected else: raise ValueError("Unknown interpolation: %s" % interpolation) + if circle: + reconstructed[~reconstruction_circle] = 0. return reconstructed * np.pi / (2 * len(th)) From 2be327815ebd6e8a6694e016871d9bcafe7ee55c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sat, 6 Jul 2013 17:28:49 +0200 Subject: [PATCH 4/8] radon: Use numpy.interp/scipy.interpolate. Scipy supports all interpolation kinds (nearest, linear) we need, while numpy supports only linear interpolation. The numpy interpolation is substantially faster, so this is used even though it complicates the code slightly. --- skimage/transform/radon_transform.py | 41 ++++++++++++++-------------- 1 file changed, 20 insertions(+), 21 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index 71548324..b99da511 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -16,6 +16,7 @@ References: from __future__ import division import numpy as np from scipy.fftpack import fftshift, fft, ifft +from scipy.interpolate import interp1d from ._warps_cy import _warp_fast from ._radon_transform import sart_projection_update from .. import util @@ -137,9 +138,10 @@ def iradon(radon_image, theta=None, output_size=None, Filter used in frequency domain filtering. Ramp filter used by default. Filters available: ramp, shepp-logan, cosine, hamming, hann Assign None to use no filter. - interpolation : str, optional (default linear) + interpolation : str, optional (default 'linear') Interpolation method used in reconstruction. - Methods available: nearest, linear. + Methods available are the same as for `scipy.interpolate.interp1d: + 'linear', 'nearest', 'zero', 'slinear', 'quadratic' and 'cubic'. circle : boolean, optional Assume the reconstructed image is zero outside the inscribed circle. Also changes the default output_size to match the behaviour of @@ -167,6 +169,10 @@ def iradon(radon_image, theta=None, output_size=None, if len(theta) != radon_image.shape[1]: raise ValueError("The given ``theta`` does not match the number of " "projections in ``radon_image``.") + interpolation_types = ('linear', 'nearest', 'zero', 'slinear', + 'quadratic', 'cubic') + if not interpolation in interpolation_types: + raise ValueError("Unknown interpolation: %s" % interpolation) if not output_size: # If output size not specified, estimate from input radon image if circle: @@ -215,7 +221,7 @@ def iradon(radon_image, theta=None, output_size=None, # Determine the center of the projections (= center of sinogram) circle_size = int(np.floor(radon_image.shape[0] / np.sqrt(2))) square_size = radon_image.shape[0] - mid_index = (square_size - circle_size) // 2 + circle_size // 2 + 1 + mid_index = (square_size - circle_size) // 2 + circle_size // 2 x = output_size y = output_size @@ -227,24 +233,17 @@ def iradon(radon_image, theta=None, output_size=None, reconstruction_circle = (xpr**2 + ypr**2) < radius**2 # Reconstruct image by interpolation - if interpolation == "nearest": - for i in range(len(theta)): - k = np.round(mid_index + ypr * np.cos(th[i]) - xpr * np.sin(th[i])) - backprojected = radon_filtered[ - ((((k > 0) & (k < n)) * k) - 1).astype(np.int), i] - reconstructed += backprojected - elif interpolation == "linear": - for i in range(len(theta)): - t = ypr * np.cos(th[i]) - xpr * np.sin(th[i]) - a = np.floor(t) - b = mid_index + a - b0 = ((((b + 1 > 0) & (b + 1 < n)) * (b + 1)) - 1).astype(np.int) - b1 = ((((b > 0) & (b < n)) * b) - 1).astype(np.int) - backprojected = (t - a) * radon_filtered[b0, i] + \ - (a - t + 1) * radon_filtered[b1, i] - reconstructed += backprojected - else: - raise ValueError("Unknown interpolation: %s" % interpolation) + for i in range(len(theta)): + t = ypr * np.cos(th[i]) - xpr * np.sin(th[i]) + x = np.arange(radon_filtered.shape[0]) - mid_index + if interpolation == 'linear': + backprojected = np.interp(t, x, radon_filtered[:, i], + left=0, right=0) + else: + interpolant = interp1d(x, radon_filtered[:, i], kind=interpolation, + bounds_error=False, fill_value=0) + backprojected = interpolant(t) + reconstructed += backprojected if circle: reconstructed[~reconstruction_circle] = 0. From b90ba783ef5bc483078a2522cd21b683a0761b69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sat, 6 Jul 2013 18:27:08 +0200 Subject: [PATCH 5/8] iradon: Only allow interpolation methods working well. Of the interpolation methods provided by scipy.interpolate.interp1d, only cubic has been found to work well with the tests in skimage (other methods are either identical to linear or nearest, or they produce bad reconstructions). --- skimage/transform/radon_transform.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index b99da511..2f2a18e7 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -139,9 +139,8 @@ def iradon(radon_image, theta=None, output_size=None, Filters available: ramp, shepp-logan, cosine, hamming, hann Assign None to use no filter. interpolation : str, optional (default 'linear') - Interpolation method used in reconstruction. - Methods available are the same as for `scipy.interpolate.interp1d: - 'linear', 'nearest', 'zero', 'slinear', 'quadratic' and 'cubic'. + Interpolation method used in reconstruction. Methods available: + 'linear', 'nearest', and 'cubic' ('cubic' is slow). circle : boolean, optional Assume the reconstructed image is zero outside the inscribed circle. Also changes the default output_size to match the behaviour of @@ -169,8 +168,7 @@ def iradon(radon_image, theta=None, output_size=None, if len(theta) != radon_image.shape[1]: raise ValueError("The given ``theta`` does not match the number of " "projections in ``radon_image``.") - interpolation_types = ('linear', 'nearest', 'zero', 'slinear', - 'quadratic', 'cubic') + interpolation_types = ('linear', 'nearest', 'cubic') if not interpolation in interpolation_types: raise ValueError("Unknown interpolation: %s" % interpolation) if not output_size: From 643a1bc1f58dc6015446d97a63977c5f4d5a4c2c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sat, 6 Jul 2013 18:29:05 +0200 Subject: [PATCH 6/8] iradon: Add test for cubic interpolation. --- skimage/transform/tests/test_radon_transform.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py index 333c298d..1d82b29a 100644 --- a/skimage/transform/tests/test_radon_transform.py +++ b/skimage/transform/tests/test_radon_transform.py @@ -124,11 +124,11 @@ def check_radon_iradon(interpolation_type, filter_type): print('\n\tmean error:', delta) if debug: _debug_plot(image, reconstructed) - if filter_type == 'ramp': - if interpolation_type == 'linear': - allowed_delta = 0.02 - else: + if filter_type in ('ramp', 'shepp-logan'): + if interpolation_type == 'nearest': allowed_delta = 0.03 + else: + allowed_delta = 0.02 else: allowed_delta = 0.05 assert delta < allowed_delta @@ -136,11 +136,12 @@ def check_radon_iradon(interpolation_type, filter_type): def test_radon_iradon(): filter_types = ["ramp", "shepp-logan", "cosine", "hamming", "hann"] - interpolation_types = ["linear", "nearest"] + interpolation_types = ['linear', 'nearest'] for interpolation_type, filter_type in \ itertools.product(interpolation_types, filter_types): yield check_radon_iradon, interpolation_type, filter_type - + # cubic interpolation is slow; only run one test for it + yield check_radon_iradon, 'cubic', 'shepp-logan' def test_iradon_angles(): """ From b3746b90902dbb47d5da2d8500973c6f00392ea9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sat, 6 Jul 2013 18:29:31 +0200 Subject: [PATCH 7/8] iradon: Cleanup by locating related code in one place. --- skimage/transform/radon_transform.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index 2f2a18e7..dfac0748 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -226,9 +226,6 @@ def iradon(radon_image, theta=None, output_size=None, [X, Y] = np.mgrid[0.0:x, 0.0:y] xpr = X - int(output_size) // 2 ypr = Y - int(output_size) // 2 - if circle: - radius = (output_size - 1) // 2 - reconstruction_circle = (xpr**2 + ypr**2) < radius**2 # Reconstruct image by interpolation for i in range(len(theta)): @@ -243,6 +240,8 @@ def iradon(radon_image, theta=None, output_size=None, backprojected = interpolant(t) reconstructed += backprojected if circle: + radius = (output_size - 1) // 2 + reconstruction_circle = (xpr**2 + ypr**2) < radius**2 reconstructed[~reconstruction_circle] = 0. return reconstructed * np.pi / (2 * len(th)) From bea50aa6082846a0673fca6c8d9f6d46668b78eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sat, 6 Jul 2013 18:55:06 +0200 Subject: [PATCH 8/8] iradon: use util.pad for sinogram padding. --- skimage/transform/radon_transform.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index dfac0748..6009c54a 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -182,16 +182,15 @@ def iradon(radon_image, theta=None, output_size=None, radon_image = _sinogram_circle_to_square(radon_image) th = (np.pi / 180.0) * theta - n = radon_image.shape[0] - img = radon_image.copy() - # resize image to next power of two for fourier analysis - # speeds up fourier and lessens artifacts - order = max(64., 2**np.ceil(np.log(2 * n) / np.log(2))) - # zero pad input image - img.resize((order, img.shape[1])) + # resize image to next power of two (but no less than 64) for + # Fourier analysis; speeds up Fourier and lessens artifacts + projection_size_padded = \ + max(64, int(2**np.ceil(np.log2(2 * radon_image.shape[0])))) + pad_width = ((0, projection_size_padded - radon_image.shape[0]), (0, 0)) + img = util.pad(radon_image, pad_width, mode='constant', constant_values=0) # Construct the Fourier filter - f = fftshift(abs(np.mgrid[-1:1:2 / order])).reshape(-1, 1) + 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 if filter == "ramp":