diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index 59e6c338..b114adad 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -72,49 +72,52 @@ def radon(image, theta=None, circle=False): 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: height, width = image.shape - diagonal = np.sqrt(height**2 + width**2) - heightpad = np.ceil(diagonal - height) - widthpad = np.ceil(diagonal - width) + 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 = 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)) - + y0 = heightpad // 2 + y1 = y0 + height + x0 = widthpad // 2 + x1 = x0 + width padded_image[y0:y1, x0:x1] = image out = np.zeros((max(padded_image.shape), len(theta))) + dh = y0 + height // 2 + dw = x0 + width // 2 - h, w = padded_image.shape - dh, dw = h // 2, w // 2 shift0 = np.array([[1, 0, -dw], [0, 1, -dh], [0, 0, 1]]) - shift1 = np.array([[1, 0, dw], [0, 1, dh], [0, 0, 1]]) def build_rotation(theta): - T = -np.deg2rad(theta) - - R = np.array([[np.cos(T), -np.sin(T), 0], - [np.sin(T), np.cos(T), 0], + T = np.deg2rad(theta) + R = np.array([[np.cos(T), np.sin(T), 0], + [-np.sin(T), np.cos(T), 0], [0, 0, 1]]) - return shift1.dot(R).dot(shift0) for i in range(len(theta)): - rotated = _warp_fast(padded_image, - np.linalg.inv(build_rotation(-theta[i]))) - - out[:, i] = rotated.sum(0)[::-1] - + rotated = _warp_fast(padded_image, build_rotation(theta[i])) + out[:, i] = rotated.sum(0) return out +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 + + def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolation="linear", circle=False): """ @@ -152,53 +155,44 @@ def iradon(radon_image, theta=None, output_size=None, Notes ----- - It applies the fourier slice theorem to reconstruct an image by + It applies the Fourier slice theorem to reconstruct an image by multiplying the frequency domain of the filter with the FFT of the projection data. This algorithm is called filtered back projection. """ if radon_image.ndim != 2: raise ValueError('The input image must be 2-D') - if theta is None: m, n = radon_image.shape theta = np.linspace(0, 180, n, endpoint=False) else: theta = np.asarray(theta) - if len(theta) != radon_image.shape[1]: raise ValueError("The given ``theta`` does not match the number of " "projections in ``radon_image``.") - - th = (np.pi / 180.0) * theta - # if output size not specified, estimate from input radon image if not output_size: + # If output size not specified, estimate from input radon image 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 + 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])) - # construct the fourier filter + # Construct the Fourier filter f = fftshift(abs(np.mgrid[-1:1:2 / order])).reshape(-1, 1) w = 2 * np.pi * f - # start from first element to avoid divide by zero + # Start from first element to avoid divide by zero if filter == "ramp": pass elif filter == "shepp-logan": @@ -213,41 +207,40 @@ def iradon(radon_image, theta=None, output_size=None, f[1:] = 1 else: raise ValueError("Unknown filter: %s" % filter) - filter_ft = np.tile(f, (1, len(theta))) - - # apply filter in fourier domain + # Apply filter in Fourier domain projection = fft(img, axis=0) * filter_ft radon_filtered = np.real(ifft(projection, axis=0)) - # resize filtered image back to original size + # Resize filtered image back to original size radon_filtered = radon_filtered[:radon_image.shape[0], :] reconstructed = np.zeros((output_size, output_size)) - mid_index = np.ceil(n / 2.0) + # 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 x = output_size y = output_size [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 + # 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])) + 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)): - t = xpr * np.sin(th[i]) - ypr * np.cos(th[i]) + 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) diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py index 79cf80c5..f7127aa5 100644 --- a/skimage/transform/tests/test_radon_transform.py +++ b/skimage/transform/tests/test_radon_transform.py @@ -2,9 +2,39 @@ from __future__ import print_function from __future__ import division import numpy as np -from numpy.testing import * +from numpy.testing import assert_raises import itertools -from skimage.transform import * +from skimage.transform import radon, iradon +from skimage.io import imread +from skimage import data_dir + + +__PHANTOM = imread(data_dir + "/phantom.png", as_grey=True)[::2, ::2] + + +def _get_phantom(): + return __PHANTOM + + +def _debug_plot(original, result, sinogram=None): + from matplotlib import pyplot as plt + imkwargs = dict(cmap='gray', interpolation='nearest') + if sinogram is None: + plt.figure(figsize=(15, 6)) + sp = 130 + else: + plt.figure(figsize=(11, 11)) + sp = 221 + plt.subplot(sp + 0) + plt.imshow(sinogram, aspect='auto', **imkwargs) + plt.subplot(sp + 1) + plt.imshow(original, **imkwargs) + plt.subplot(sp + 2) + plt.imshow(result, vmin=original.min(), vmax=original.max(), **imkwargs) + plt.subplot(sp + 3) + plt.imshow(result - original, **imkwargs) + plt.colorbar() + plt.show() def rescale(x): @@ -14,35 +44,100 @@ def rescale(x): return x -def test_radon_iradon(): - size = 100 +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) + image[(shape[0] // 2, shape[1] // 2)] = 1. + # Calculate the sinogram + theta = np.linspace(0., 180., max(shape), endpoint=False) + sinogram = radon(image, theta=theta, circle=circle) + # The sinogram should be a straight, horizontal line + sinogram_max = np.argmax(sinogram, axis=0) + print(sinogram_max) + assert np.std(sinogram_max) < 1e-6 + + +def test_radon_center(): + shapes = [(16, 16), (17, 17)] + circles = [False, True] + for shape, circle in itertools.product(shapes, circles): + yield check_radon_center, shape, circle + rectangular_shapes = [(32, 16), (33, 17)] + for shape in rectangular_shapes: + yield check_radon_center, shape, False + + +def check_iradon_center(size, theta, circle): debug = False - image = np.tri(size) + np.tri(size)[::-1] - for filter_type in ["ramp", "shepp-logan", "cosine", "hamming", "hann"]: - reconstructed = iradon(radon(image), filter=filter_type) + # Create a test sinogram corresponding to a single projection + # with a single non-zero pixel at the rotation center + if circle: + sinogram = np.zeros((size, 1), dtype=np.float) + sinogram[size // 2, 0] = 1. + else: + diagonal = int(np.ceil(np.sqrt(2) * size)) + sinogram = np.zeros((diagonal, 1), dtype=np.float) + sinogram[sinogram.shape[0] // 2, 0] = 1. + maxpoint = np.unravel_index(np.argmax(sinogram), sinogram.shape) + print('shape of generated sinogram', sinogram.shape) + print('maximum in generated sinogram', maxpoint) + # Compare reconstructions for theta=angle and theta=angle + 180; + # these should be exactly equal + reconstruction = iradon(sinogram, theta=[theta], circle=circle) + reconstruction_opposite = iradon(sinogram, theta=[theta + 180], + circle=circle) + print('rms deviance:', + np.sqrt(np.mean((reconstruction_opposite - reconstruction)**2))) + if debug: + import matplotlib.pyplot as plt + imkwargs = dict(cmap='gray', interpolation='nearest') + plt.figure() + plt.subplot(221) + plt.imshow(sinogram, **imkwargs) + plt.subplot(222) + plt.imshow(reconstruction_opposite - reconstruction, **imkwargs) + plt.subplot(223) + plt.imshow(reconstruction, **imkwargs) + plt.subplot(224) + plt.imshow(reconstruction_opposite, **imkwargs) + plt.show() - image = rescale(image) - reconstructed = rescale(reconstructed) - delta = np.mean(np.abs(image - reconstructed)) + assert np.allclose(reconstruction, reconstruction_opposite) - if debug: - print(delta) - import matplotlib.pyplot as plt - f, (ax1, ax2) = plt.subplots(1, 2) - ax1.imshow(image, cmap=plt.cm.gray) - ax2.imshow(reconstructed, cmap=plt.cm.gray) - plt.show() - assert delta < 0.05 +def test_iradon_center(): + sizes = [16, 17] + thetas = [0, 90] + circles = [False, True] + for size, theta, circle in itertools.product(sizes, thetas, circles): + yield check_iradon_center, size, theta, circle - reconstructed = iradon(radon(image), filter="ramp", - interpolation="nearest") - delta = np.mean(abs(image - reconstructed)) - assert delta < 0.05 - size = 20 - image = np.tri(size) + np.tri(size)[::-1] - reconstructed = iradon(radon(image), filter="ramp", - interpolation="nearest") + +def check_radon_iradon(interpolation_type, filter_type): + debug = False + image = _get_phantom() + reconstructed = iradon(radon(image), filter=filter_type, + interpolation=interpolation_type) + delta = np.mean(np.abs(image - reconstructed)) + print('\n\tmean error:', delta) + if debug: + _debug_plot(image, reconstructed) + if filter_type == 'ramp': + if interpolation_type == 'linear': + allowed_delta = 0.02 + else: + allowed_delta = 0.03 + else: + allowed_delta = 0.05 + assert delta < allowed_delta + + +def test_radon_iradon(): + filter_types = ["ramp", "shepp-logan", "cosine", "hamming", "hann"] + interpolation_types = ["linear", "nearest"] + for interpolation_type, filter_type in \ + itertools.product(interpolation_types, filter_types): + yield check_radon_iradon, interpolation_type, filter_type def test_iradon_angles(): @@ -73,32 +168,29 @@ def test_iradon_angles(): assert delta_80 > delta_200 -def test_radon_minimal(): - """ - Test for small images for various angles - """ - thetas = [np.arange(180)] - for theta in thetas: - a = np.zeros((3, 3)) - a[1, 1] = 1 - p = radon(a, theta) - reconstructed = iradon(p, theta) - reconstructed /= np.max(reconstructed) - assert np.all(abs(a - reconstructed) < 0.4) +def check_radon_iradon_minimal(shape, slices): + debug = False + theta = np.arange(180) + image = np.zeros(shape, dtype=np.float) + image[slices] = 1. + sinogram = radon(image, theta) + reconstructed = iradon(sinogram, theta) + print('\n\tMaximum deviation:', np.max(np.abs(image - reconstructed))) + if debug: + _debug_plot(image, reconstructed, sinogram) + if image.sum() == 1: + assert (np.unravel_index(np.argmax(reconstructed), image.shape) + == np.unravel_index(np.argmax(image), image.shape)) - b = np.zeros((4, 4)) - b[1:3, 1:3] = 1 - p = radon(b, theta) - reconstructed = iradon(p, theta) - reconstructed /= np.max(reconstructed) - assert np.all(abs(b - reconstructed) < 0.4) - c = np.zeros((5, 5)) - c[1:3, 1:3] = 1 - p = radon(c, theta) - reconstructed = iradon(p, theta) - reconstructed /= np.max(reconstructed) - assert np.all(abs(c - reconstructed) < 0.4) +def test_radon_iradon_minimal(): + shapes = [(3, 3), (4, 4), (5, 5)] + for shape in shapes: + c0, c1 = shape[0] // 2, shape[1] // 2 + coordinates = itertools.product((c0 - 1, c0, c0 + 1), + (c1 - 1, c1, c1 + 1)) + for coordinate in coordinates: + yield check_radon_iradon_minimal, shape, coordinate def test_reconstruct_with_wrong_angles(): @@ -144,40 +236,64 @@ def test_radon_circle(): assert np.all(relative_error < 3e-3) +def check_sinogram_circle_to_square(size): + from skimage.transform.radon_transform import _sinogram_circle_to_square + image = _random_circle((size, size)) + theta = np.linspace(0., 180., size, False) + sinogram_circle = radon(image, theta, circle=True) + argmax_shape = lambda a: np.unravel_index(np.argmax(a), a.shape) + print('\n\targmax of circle:', argmax_shape(sinogram_circle)) + sinogram_square = radon(image, theta, circle=False) + print('\targmax of square:', argmax_shape(sinogram_square)) + sinogram_circle_to_square = _sinogram_circle_to_square(sinogram_circle) + print('\targmax of circle to square:', + argmax_shape(sinogram_circle_to_square)) + error = abs(sinogram_square - sinogram_circle_to_square) + print(np.mean(error), np.max(error)) + assert (argmax_shape(sinogram_square) + == argmax_shape(sinogram_circle_to_square)) + + +def test_sinogram_circle_to_square(): + for size in (50, 51): + yield check_sinogram_circle_to_square, size + + +def check_radon_iradon_circle(interpolation, shape, output_size): + # Forward and inverse radon on synthetic data + image = _random_circle(shape) + radius = min(shape) // 2 + 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) + + def test_radon_iradon_circle(): shape = (61, 79) - radius = min(shape) // 2 - image = _random_circle(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) - # 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) + yield check_radon_iradon_circle, interpolation, shape, output_size if __name__ == "__main__":