From 2bfaafd9e2bdc5470c1c38a1cf568f49d0f00077 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sat, 13 Jul 2013 19:31:43 +0200 Subject: [PATCH 1/8] Radon transform: Redefine projection center for sinograms. This definition is chosen because it is simple to express in the documentation. No changes in accuracy are to be expected, but comparing sinograms and reconstructions before and after this commit will give different results in the cases where ``circle=False`` for ``radon`` or ``iradon``. --- skimage/transform/radon_transform.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index 1ebd26c6..9fe1c6e1 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -82,12 +82,15 @@ def radon(image, theta=None, circle=False): else: diagonal = np.sqrt(2) * max(image.shape) pad = [int(np.ceil(diagonal - s)) for s in image.shape] - pad_width = [(p // 2, p - p // 2) for p in pad] + new_center = [(s + p) // 2 for s, p in zip(image.shape, pad)] + old_center = [s // 2 for s in image.shape] + pad_before = [nc - oc for oc, nc in zip(old_center, new_center)] + pad_width = [(pb, p - pb) for pb, p in zip(pad_before, pad)] padded_image = util.pad(image, pad_width, mode='constant', constant_values=0) out = np.zeros((max(padded_image.shape), len(theta))) - dh = pad[0] // 2 + image.shape[0] // 2 - dw = pad[1] // 2 + image.shape[1] // 2 + dh = padded_image.shape[0] // 2 + dw = padded_image.shape[1] // 2 shift0 = np.array([[1, 0, -dw], [0, 1, -dh], @@ -112,7 +115,10 @@ def radon(image, theta=None, circle=False): def _sinogram_circle_to_square(sinogram): diagonal = int(np.ceil(np.sqrt(2) * sinogram.shape[0])) pad = diagonal - sinogram.shape[0] - pad_width = ((pad // 2, pad - pad // 2), (0, 0)) + old_center = sinogram.shape[0] // 2 + new_center = diagonal // 2 + pad_before = new_center - old_center + pad_width = ((pad_before, pad - pad_before), (0, 0)) return util.pad(sinogram, pad_width, mode='constant', constant_values=0) @@ -216,9 +222,7 @@ def iradon(radon_image, theta=None, output_size=None, radon_filtered = radon_filtered[:radon_image.shape[0], :] reconstructed = np.zeros((output_size, output_size)) # 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 + mid_index = radon_image.shape[0] // 2 x = output_size y = output_size From 47b6d0c5a625aae4fda3348944c6738e3331588f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sat, 13 Jul 2013 19:34:22 +0200 Subject: [PATCH 2/8] Radon transform: Document rotation axis location. --- skimage/transform/radon_transform.py | 29 ++++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index 9fe1c6e1..37505bfa 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -33,7 +33,8 @@ def radon(image, theta=None, circle=False): Parameters ---------- image : array_like, dtype=float - Input image. + Input image. The rotation axis will be located in the pixel with + indices ``(image.shape[0] // 2, image.shape[1] // 2)``. theta : array_like, dtype=float, optional (default np.arange(180)) Projection angles (in degrees). circle : boolean, optional @@ -43,8 +44,10 @@ def radon(image, theta=None, circle=False): Returns ------- - output : ndarray - Radon transform (sinogram). + radon_image : ndarray + Radon transform (sinogram). The tomography rotation axis will lie + at the pixel index ``radon_image.shape[0] // 2`` along the 0th + dimension of ``radon_image``. Raises ------ @@ -134,7 +137,10 @@ def iradon(radon_image, theta=None, output_size=None, ---------- radon_image : array_like, dtype=float Image containing radon transform (sinogram). Each column of - the image corresponds to a projection along a different angle. + the image corresponds to a projection along a different angle. The + tomography rotation axis should lie at the pixel index + ``radon_image.shape[0] // 2`` along the 0th dimension of + ``radon_image``. theta : array_like, dtype=float, optional Reconstruction angles (in degrees). Default: m angles evenly spaced between 0 and 180 (if the shape of `radon_image` is (N, M)). @@ -155,7 +161,9 @@ def iradon(radon_image, theta=None, output_size=None, Returns ------- output : ndarray - Reconstructed image. + Reconstructed image. The rotation axis will be located in the pixel + with indices + ``(reconstructed.shape[0] // 2, reconstructed.shape[1] // 2)``. Notes ----- @@ -318,7 +326,10 @@ def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None, ---------- radon_image : 2D array, dtype=float Image containing radon transform (sinogram). Each column of - the image corresponds to a projection along a different angle. + the image corresponds to a projection along a different angle. The + tomography rotation axis should lie at the pixel index + ``radon_image.shape[0] // 2`` along the 0th dimension of + ``radon_image``. theta : 1D array, dtype=float, optional Reconstruction angles (in degrees). Default: m angles evenly spaced between 0 and 180 (if the shape of `radon_image` is (N, M)). @@ -340,8 +351,10 @@ def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None, Returns ------- - output : ndarray - Reconstructed image. + reconstructed : ndarray + Reconstructed image. The rotation axis will be located in the pixel + with indices + ``(reconstructed.shape[0] // 2, reconstructed.shape[1] // 2)``. Notes ----- From 41c6f6d740a5c8936d3319bcf13dd1f204d4fd01 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sat, 13 Jul 2013 19:47:38 +0200 Subject: [PATCH 3/8] radon: Reduce duplication; simplifications. --- skimage/transform/radon_transform.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index 37505bfa..da01d1b4 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -68,6 +68,7 @@ def radon(image, theta=None, circle=False): if not np.all(reconstruction_circle | (image == 0)): raise ValueError('Image must be zero outside the reconstruction' ' circle') + # Crop image to make it square slices = [] for d in (0, 1): if image.shape[d] > min(image.shape): @@ -79,9 +80,6 @@ def radon(image, theta=None, circle=False): slices.append(slice(None)) slices = tuple(slices) padded_image = image[slices] - out = np.zeros((min(padded_image.shape), len(theta))) - dh = padded_image.shape[0] // 2 - dw = padded_image.shape[1] // 2 else: diagonal = np.sqrt(2) * max(image.shape) pad = [int(np.ceil(diagonal - s)) for s in image.shape] @@ -91,15 +89,16 @@ def radon(image, theta=None, circle=False): pad_width = [(pb, p - pb) for pb, p in zip(pad_before, pad)] padded_image = util.pad(image, pad_width, mode='constant', constant_values=0) - out = np.zeros((max(padded_image.shape), len(theta))) - dh = padded_image.shape[0] // 2 - dw = padded_image.shape[1] // 2 + # padded_image is always square + assert padded_image.shape[0] == padded_image.shape[1] + radon_image = np.zeros((padded_image.shape[0], len(theta))) + center = padded_image.shape[0] // 2 - shift0 = np.array([[1, 0, -dw], - [0, 1, -dh], + shift0 = np.array([[1, 0, -center], + [0, 1, -center], [0, 0, 1]]) - shift1 = np.array([[1, 0, dw], - [0, 1, dh], + shift1 = np.array([[1, 0, center], + [0, 1, center], [0, 0, 1]]) def build_rotation(theta): @@ -111,8 +110,8 @@ def radon(image, theta=None, circle=False): for i in range(len(theta)): rotated = _warp_fast(padded_image, build_rotation(theta[i])) - out[:, i] = rotated.sum(0) - return out + radon_image[:, i] = rotated.sum(0) + return radon_image def _sinogram_circle_to_square(sinogram): @@ -160,7 +159,7 @@ def iradon(radon_image, theta=None, output_size=None, Returns ------- - output : ndarray + reconstructed : ndarray Reconstructed image. The rotation axis will be located in the pixel with indices ``(reconstructed.shape[0] // 2, reconstructed.shape[1] // 2)``. From 288ee694834ee20c4b750085a40ef598fe807e2d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sat, 13 Jul 2013 19:50:20 +0200 Subject: [PATCH 4/8] radon: Correct docstring of order_angles_golden_ratio. --- skimage/transform/radon_transform.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index da01d1b4..6ae6e1b9 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -269,9 +269,11 @@ def order_angles_golden_ratio(theta): Returns ------- - indices : 1D array of unsigned integers - Indices into ``theta`` such that ``theta[indices]`` gives the - approximate golden ratio ordering of the projections. + indices_generator : generator yielding unsigned integers + The returned generator yields indices into ``theta`` such that + ``theta[indices]`` gives the approximate golden ratio ordering + of the projections. In total, ``len(theta)`` indices are yielded. + All non-negative integers < ``len(theta)`` are yielded exactly once. Notes ----- From 462173a53adb487b1200ed9dcc2639d7cee32fc1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sat, 13 Jul 2013 20:06:59 +0200 Subject: [PATCH 5/8] Radon transform: Include boundary in reconstruction circle. A test criterion needed to be relaxed slightly to have tests still passing. This is ok, as the reconstruction circle is now larger, meaning larger errors should be expected. Moreover, the test in question uses random data, and changing the seed causes greater changes in accuracy than the amount the test criterion was relaxed by. --- skimage/transform/radon_transform.py | 10 ++++------ skimage/transform/tests/test_radon_transform.py | 6 +++--- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index 6ae6e1b9..57935d43 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -64,7 +64,7 @@ def radon(image, theta=None, circle=False): 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 + + (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') @@ -231,9 +231,7 @@ def iradon(radon_image, theta=None, output_size=None, # Determine the center of the projections (= center of sinogram) mid_index = radon_image.shape[0] // 2 - x = output_size - y = output_size - [X, Y] = np.mgrid[0.0:x, 0.0:y] + [X, Y] = np.mgrid[0:output_size, 0:output_size] xpr = X - int(output_size) // 2 ypr = Y - int(output_size) // 2 @@ -250,8 +248,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 + radius = output_size // 2 + reconstruction_circle = (xpr**2 + ypr**2) <= radius**2 reconstructed[~reconstruction_circle] = 0. return reconstructed * np.pi / (2 * len(th)) diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py index 1d82b29a..32d8780c 100644 --- a/skimage/transform/tests/test_radon_transform.py +++ b/skimage/transform/tests/test_radon_transform.py @@ -210,7 +210,7 @@ def _random_circle(shape): 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. + image[r > radius] = 0. return image @@ -236,7 +236,7 @@ def test_radon_circle(): 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) + assert np.all(relative_error < 3.2e-3) def check_sinogram_circle_to_square(size): @@ -284,7 +284,7 @@ def check_radon_iradon_circle(interpolation, shape, output_size): # 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. + reconstruction_rectangle[r > radius] = 0. print(reconstruction_circle.shape) print(reconstruction_rectangle.shape) np.allclose(reconstruction_rectangle, reconstruction_circle) From 60444ee4d333c4079e489abbfd58b4064eb9d352 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sat, 13 Jul 2013 20:11:49 +0200 Subject: [PATCH 6/8] Radon transform: PEP8 fixes. --- skimage/transform/radon_transform.py | 6 +++--- skimage/transform/tests/test_radon_transform.py | 1 + 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index 57935d43..db43090b 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -374,9 +374,9 @@ def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None, ---------- .. [1] AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging", IEEE Press 1988. - .. [2] AH Andersen, AC Kak, "Simultaneous algebraic reconstruction technique - (SART): a superior implementation of the ART algorithm", Ultrasonic - Imaging 6 pp 81--94 (1984) + .. [2] AH Andersen, AC Kak, "Simultaneous algebraic reconstruction + technique (SART): a superior implementation of the ART algorithm", + Ultrasonic Imaging 6 pp 81--94 (1984) .. [3] S Kaczmarz, "Angenäherte auflösung von systemen linearer gleichungen", Bulletin International de l’Academie Polonaise des Sciences et des Lettres 35 pp 355--357 (1937) diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py index 32d8780c..6fbc3f62 100644 --- a/skimage/transform/tests/test_radon_transform.py +++ b/skimage/transform/tests/test_radon_transform.py @@ -143,6 +143,7 @@ def test_radon_iradon(): # cubic interpolation is slow; only run one test for it yield check_radon_iradon, 'cubic', 'shepp-logan' + def test_iradon_angles(): """ Test with different number of projections 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 7/8] 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 From d5b72f91ab4c5db34aa29c2f44d730e2cd0f3e4d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sun, 14 Jul 2013 17:40:17 +0200 Subject: [PATCH 8/8] iradon: Do not suppress 0 frequency for filter=None. --- skimage/transform/radon_transform.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index fc827c8b..101d09c6 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -218,7 +218,7 @@ def iradon(radon_image, theta=None, output_size=None, elif filter == "hann": fourier_filter *= (1 + np.cos(omega / 2)) / 2 elif filter is None: - f[1:] = 1 + fourier_filter[:] = 1 else: raise ValueError("Unknown filter: %s" % filter) # Apply filter in Fourier domain