From b87264844199c3ba3d5b50d410799ea03038a54d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sun, 26 May 2013 14:00:54 +0200 Subject: [PATCH 1/8] Add reconstruction circle option to transform.radon. --- skimage/transform/radon_transform.py | 56 +++++++++++++++++++++------- 1 file changed, 43 insertions(+), 13 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index 0213b2bc..ab312158 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -20,7 +20,7 @@ from ._warps_cy import _warp_fast __all__ = ["radon", "iradon"] -def radon(image, theta=None): +def radon(image, theta=None, circle=False): """ Calculates the radon transform of an image given specified projection angles. @@ -31,31 +31,61 @@ def radon(image, theta=None): Input image. theta : array_like, dtype=float, optional (default np.arange(180)) Projection angles (in degrees). + circle : boolean, optional (default False) + Assume image is zero outside the inscribed circle, making the + width of each projection (the first dimension of the sinogram) + equal to min(image.shape). Returns ------- output : ndarray Radon transform (sinogram). + Raises + ------ + ValueError + If called with `circle=True` and image != 0 outside the inscribed + circle """ if image.ndim != 2: raise ValueError('The input image must be 2-D') if theta is None: theta = np.arange(180) - height, width = image.shape - diagonal = np.sqrt(height**2 + width**2) - heightpad = np.ceil(diagonal - height) - widthpad = np.ceil(diagonal - width) - padded_image = np.zeros((int(height + heightpad), - int(width + widthpad))) - y0, y1 = int(np.ceil(heightpad / 2)), \ - int((np.ceil(heightpad / 2) + height)) - x0, x1 = int((np.ceil(widthpad / 2))), \ - int((np.ceil(widthpad / 2) + width)) + if circle: + radius = min(image.shape) // 2 + c0, c1 = np.ogrid[0:image.shape[0], 0:image.shape[1]] + reconstruction_circle = ((c0 - image.shape[0] // 2)**2 + + (c1 - image.shape[1] // 2)**2) < radius**2 + if not np.all(reconstruction_circle | (image == 0)): + raise ValueError('image must be zero outside the reconstruction' + + ' circle') + slices = [] + for d in (0, 1): + if image.shape[d] > min(image.shape): + excess = image.shape[d] - min(image.shape) + slices.append(slice(int(np.ceil(excess / 2)), + int(np.ceil(excess / 2) + + min(image.shape)))) + else: + slices.append(slice(None)) + slices = tuple(slices) + padded_image = image[slices] + out = np.zeros((min(padded_image.shape), len(theta))) + else: + height, width = image.shape + diagonal = np.sqrt(height**2 + width**2) + heightpad = np.ceil(diagonal - height) + widthpad = np.ceil(diagonal - width) + padded_image = np.zeros((int(height + heightpad), + int(width + widthpad))) + y0, y1 = int(np.ceil(heightpad / 2)), \ + int((np.ceil(heightpad / 2) + height)) + x0, x1 = int((np.ceil(widthpad / 2))), \ + int((np.ceil(widthpad / 2) + width)) - padded_image[y0:y1, x0:x1] = image - out = np.zeros((max(padded_image.shape), len(theta))) + padded_image[y0:y1, x0:x1] = image + out = np.zeros((max(padded_image.shape), len(theta))) h, w = padded_image.shape dh, dw = h // 2, w // 2 From aa8f1b7c98cf49652e9edac4a86fc0c0941f574d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sun, 26 May 2013 14:03:04 +0200 Subject: [PATCH 2/8] Tests for reconstruction circle option for transform.radon. --- .../transform/tests/test_radon_transform.py | 27 +++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py index 3b2f19dc..2e1a016f 100644 --- a/skimage/transform/tests/test_radon_transform.py +++ b/skimage/transform/tests/test_radon_transform.py @@ -1,4 +1,5 @@ from __future__ import print_function +from __future__ import division import numpy as np from numpy.testing import * @@ -103,6 +104,32 @@ def test_reconstruct_with_wrong_angles(): iradon(p, theta=[0, 1, 2]) assert_raises(ValueError, iradon, p, theta=[0, 1, 2, 3]) +def test_radon_circle(): + a = np.ones((10, 10)) + assert_raises(ValueError, radon, a, circle=True) + + # Synthetic data, circular symmetry + shape = (61, 79) + c0, c1 = np.ogrid[0:shape[0], 0:shape[1]] + r = np.sqrt((c0 - shape[0] // 2)**2 + (c1 - shape[1] // 2)**2) + radius = min(shape) // 2 + image = np.clip(radius - r, 0, np.inf) + image = rescale(image) + angles = np.linspace(0, 180, min(shape), endpoint=False) + sinogram = radon(image, theta=angles, circle=True) + assert np.all(sinogram.std(axis=1) < 1e-2) + + # Synthetic data, random + np.random.seed(98312871) + image = np.random.rand(*shape) + image[r >= radius] = 0. + sinogram = radon(image, theta=angles, circle=True) + mass = sinogram.sum(axis=0) + average_mass = mass.mean() + relative_error = np.abs(mass - average_mass) / average_mass + print(relative_error.max(), relative_error.mean()) + assert np.all(relative_error < 3e-3) + if __name__ == "__main__": run_module_suite() From 5876d80414c8bd059b5817a51da2d6d2962a236d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sun, 26 May 2013 19:26:05 +0200 Subject: [PATCH 3/8] Add reconstruction circle option to transform.iradon. --- skimage/transform/radon_transform.py | 35 ++++++++++++++++++++++++---- 1 file changed, 30 insertions(+), 5 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index ab312158..e44940ec 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -116,7 +116,7 @@ def radon(image, theta=None, circle=False): def iradon(radon_image, theta=None, output_size=None, - filter="ramp", interpolation="linear"): + filter="ramp", interpolation="linear", circle=False): """ Inverse radon transform. @@ -140,6 +140,10 @@ def iradon(radon_image, theta=None, output_size=None, interpolation : str, optional (default linear) Interpolation method used in reconstruction. Methods available: nearest, linear. + circle : boolean, optional (default False) + Assume the reconstructed image is zero outside the inscribed circle. + Also changes the default output_size to match the behaviour of + ``radon`` called with circle=True. Returns ------- @@ -169,7 +173,19 @@ def iradon(radon_image, theta=None, output_size=None, th = (np.pi / 180.0) * theta # if output size not specified, estimate from input radon image if not output_size: - output_size = int(np.floor(np.sqrt((radon_image.shape[0])**2 / 2.0))) + if circle: + output_size = radon_image.shape[0] + else: + output_size = int(np.floor(np.sqrt((radon_image.shape[0])**2 + / 2.0))) + if circle: + radon_size = int(np.ceil(np.sqrt(2) * radon_image.shape[0])) + radon_image_padded = np.zeros((radon_size, radon_image.shape[1])) + radon_pad = (radon_size - radon_image.shape[0]) // 2 + radon_image_padded[radon_pad:radon_pad + radon_image.shape[0], :] \ + = radon_image + radon_image = radon_image_padded + n = radon_image.shape[0] img = radon_image.copy() @@ -215,12 +231,19 @@ def iradon(radon_image, theta=None, output_size=None, 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 if interpolation == "nearest": for i in range(len(theta)): k = np.round(mid_index + xpr * np.sin(th[i]) - ypr * np.cos(th[i])) - reconstructed += radon_filtered[ + 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)): @@ -229,9 +252,11 @@ def iradon(radon_image, theta=None, output_size=None, 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) - reconstructed += (t - a) * radon_filtered[b0, i] + \ + 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) From 2118f1d97db620e0251f55a3f5e2c304d3e156e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sun, 26 May 2013 19:26:39 +0200 Subject: [PATCH 4/8] Tests for reconstruction circle mode in transform.iradon. --- .../transform/tests/test_radon_transform.py | 40 +++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py index 2e1a016f..e6faa78f 100644 --- a/skimage/transform/tests/test_radon_transform.py +++ b/skimage/transform/tests/test_radon_transform.py @@ -3,6 +3,7 @@ from __future__ import division import numpy as np from numpy.testing import * +import itertools from skimage.transform import * @@ -130,6 +131,45 @@ def test_radon_circle(): print(relative_error.max(), relative_error.mean()) assert np.all(relative_error < 3e-3) +def test_radon_iradon_circle(): + shape = (61, 79) + # Synthetic random data, zero outside reconstruction circle + image = np.random.rand(*shape) + interpolations = ('nearest', 'linear') + output_sizes = (None, min(shape), max(shape), 97) + + for interpolation, output_size in itertools.product(interpolations, + output_sizes): + print('interpolation =', interpolation) + print('output_size =', output_size) + c0, c1 = np.ogrid[0:shape[0], 0:shape[1]] + r = np.sqrt((c0 - shape[0] // 2)**2 + (c1 - shape[1] // 2)**2) + radius = min(shape) // 2 + image[r >= radius] = 0. + # Forward and inverse radon on synthetic data + sinogram_rectangle = radon(image, circle=False) + reconstruction_rectangle = iradon(sinogram_rectangle, + output_size=output_size, + interpolation=interpolation, + circle=False) + sinogram_circle = radon(image, circle=True) + reconstruction_circle = iradon(sinogram_circle, + output_size=output_size, + interpolation=interpolation, + circle=True) + # Crop rectangular reconstruction to match circle=True reconstruction + width = reconstruction_circle.shape[0] + excess = int(np.ceil((reconstruction_rectangle.shape[0] - width) / 2)) + s = np.s_[excess:width + excess, excess:width + excess] + reconstruction_rectangle = reconstruction_rectangle[s] + # Find the reconstruction circle, set reconstruction to zero outside + c0, c1 = np.ogrid[0:width, 0:width] + r = np.sqrt((c0 - width // 2)**2 + (c1 - width // 2)**2) + reconstruction_rectangle[r >= radius] = 0. + print(reconstruction_circle.shape) + print(reconstruction_rectangle.shape) + np.allclose(reconstruction_rectangle, reconstruction_circle) + if __name__ == "__main__": run_module_suite() From 8208293c93c67c9cf620010070160bfc7c6f1ffd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sun, 26 May 2013 21:15:23 +0200 Subject: [PATCH 5/8] Correct documentation and string formatting. --- skimage/transform/radon_transform.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index e44940ec..9ce1f28e 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -34,7 +34,7 @@ def radon(image, theta=None, circle=False): circle : boolean, optional (default False) Assume image is zero outside the inscribed circle, making the width of each projection (the first dimension of the sinogram) - equal to min(image.shape). + equal to ``min(image.shape)``. Returns ------- @@ -44,7 +44,7 @@ def radon(image, theta=None, circle=False): Raises ------ ValueError - If called with `circle=True` and image != 0 outside the inscribed + If called with ``circle=True`` and ``image != 0`` outside the inscribed circle """ if image.ndim != 2: @@ -58,8 +58,8 @@ def radon(image, theta=None, circle=False): reconstruction_circle = ((c0 - image.shape[0] // 2)**2 + (c1 - image.shape[1] // 2)**2) < radius**2 if not np.all(reconstruction_circle | (image == 0)): - raise ValueError('image must be zero outside the reconstruction' - + ' circle') + raise ValueError('Image must be zero outside the reconstruction' + ' circle') slices = [] for d in (0, 1): if image.shape[d] > min(image.shape): @@ -143,7 +143,7 @@ def iradon(radon_image, theta=None, output_size=None, circle : boolean, optional (default False) Assume the reconstructed image is zero outside the inscribed circle. Also changes the default output_size to match the behaviour of - ``radon`` called with circle=True. + ``radon`` called with ``circle=True``. Returns ------- From c9ae913d325fef9f86d3f0571efe316fab56b042 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Mon, 27 May 2013 20:42:59 +0200 Subject: [PATCH 6/8] PEP8 style changes for transform.radon_transform and its tests --- skimage/transform/radon_transform.py | 12 ++++++------ skimage/transform/tests/test_radon_transform.py | 2 ++ 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index 9ce1f28e..f31a0263 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -79,10 +79,10 @@ def radon(image, theta=None, circle=False): widthpad = np.ceil(diagonal - width) padded_image = np.zeros((int(height + heightpad), int(width + widthpad))) - y0, y1 = int(np.ceil(heightpad / 2)), \ - int((np.ceil(heightpad / 2) + height)) - x0, x1 = int((np.ceil(widthpad / 2))), \ - int((np.ceil(widthpad / 2) + width)) + y0 = int(np.ceil(heightpad / 2)) + y1 = int((np.ceil(heightpad / 2) + height)) + x0 = int((np.ceil(widthpad / 2))) + x1 = int((np.ceil(widthpad / 2) + width)) padded_image[y0:y1, x0:x1] = image out = np.zeros((max(padded_image.shape), len(theta))) @@ -209,7 +209,7 @@ def iradon(radon_image, theta=None, output_size=None, f[1:] = f[1:] * (0.54 + 0.46 * np.cos(w[1:])) elif filter == "hann": f[1:] = f[1:] * (1 + np.cos(w[1:])) / 2 - elif filter == None: + elif filter is None: f[1:] = 1 else: raise ValueError("Unknown filter: %s" % filter) @@ -253,7 +253,7 @@ def iradon(radon_image, theta=None, output_size=None, 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] + (a - t + 1) * radon_filtered[b1, i] if circle: backprojected[~reconstruction_circle] = 0. reconstructed += backprojected diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py index e6faa78f..9b9dc2b0 100644 --- a/skimage/transform/tests/test_radon_transform.py +++ b/skimage/transform/tests/test_radon_transform.py @@ -105,6 +105,7 @@ def test_reconstruct_with_wrong_angles(): iradon(p, theta=[0, 1, 2]) assert_raises(ValueError, iradon, p, theta=[0, 1, 2, 3]) + def test_radon_circle(): a = np.ones((10, 10)) assert_raises(ValueError, radon, a, circle=True) @@ -131,6 +132,7 @@ def test_radon_circle(): print(relative_error.max(), relative_error.mean()) assert np.all(relative_error < 3e-3) + def test_radon_iradon_circle(): shape = (61, 79) # Synthetic random data, zero outside reconstruction circle From bec5a20441a945d80ed92968756f2715faab44f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Mon, 27 May 2013 20:43:51 +0200 Subject: [PATCH 7/8] Style improvements in docstrings in transform.radon_transform --- skimage/transform/radon_transform.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index f31a0263..59e6c338 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -31,7 +31,7 @@ def radon(image, theta=None, circle=False): Input image. theta : array_like, dtype=float, optional (default np.arange(180)) Projection angles (in degrees). - circle : boolean, optional (default False) + circle : boolean, optional Assume image is zero outside the inscribed circle, making the width of each projection (the first dimension of the sinogram) equal to ``min(image.shape)``. @@ -140,7 +140,7 @@ def iradon(radon_image, theta=None, output_size=None, interpolation : str, optional (default linear) Interpolation method used in reconstruction. Methods available: nearest, linear. - circle : boolean, optional (default False) + circle : boolean, optional Assume the reconstructed image is zero outside the inscribed circle. Also changes the default output_size to match the behaviour of ``radon`` called with ``circle=True``. From d5629b12656ef63bb6a3cd0c38f68ccde751eb47 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Mon, 27 May 2013 21:47:56 +0200 Subject: [PATCH 8/8] Update contributors --- CONTRIBUTORS.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 1d018cb0..f395b68a 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -139,3 +139,6 @@ - Xavier Moles Lopez Color separation (color deconvolution) for several stainings. + +- Jostein Bø Fløystad + Reconstruction circle mode for Radon transform