494 lines
19 KiB
Python
494 lines
19 KiB
Python
import itertools
|
|
import pytest
|
|
|
|
import numpy as np
|
|
from skimage.data import shepp_logan_phantom
|
|
from skimage.transform import radon, iradon, iradon_sart, rescale
|
|
|
|
from skimage._shared.utils import convert_to_float
|
|
from skimage._shared import testing
|
|
from skimage._shared.testing import test_parallel
|
|
from skimage._shared._warnings import expected_warnings
|
|
|
|
|
|
PHANTOM = shepp_logan_phantom()[::2, ::2]
|
|
PHANTOM = rescale(PHANTOM, 0.5, order=1,
|
|
mode='constant', anti_aliasing=False, multichannel=False)
|
|
|
|
|
|
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_intensity(x):
|
|
x = x.astype(float)
|
|
x -= x.min()
|
|
x /= x.max()
|
|
return x
|
|
|
|
|
|
def test_iradon_bias_circular_phantom():
|
|
"""
|
|
test that a uniform circular phantom has a small reconstruction bias
|
|
"""
|
|
pixels = 128
|
|
xy = np.arange(-pixels / 2, pixels / 2) + 0.5
|
|
x, y = np.meshgrid(xy, xy)
|
|
image = x**2 + y**2 <= (pixels/4)**2
|
|
|
|
theta = np.linspace(0., 180., max(image.shape), endpoint=False)
|
|
sinogram = radon(image, theta=theta)
|
|
|
|
reconstruction_fbp = iradon(sinogram, theta=theta)
|
|
error = reconstruction_fbp - image
|
|
|
|
tol = 5e-5
|
|
roi_err = np.abs(np.mean(error))
|
|
assert roi_err < tol
|
|
|
|
|
|
def check_radon_center(shape, circle, dtype, preserve_range):
|
|
# Create a test image with only a single non-zero pixel at the origin
|
|
image = np.zeros(shape, dtype=dtype)
|
|
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,
|
|
preserve_range=preserve_range)
|
|
# 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
|
|
|
|
|
|
@testing.parametrize("shape", [(16, 16), (17, 17)])
|
|
@testing.parametrize("circle", [False, True])
|
|
@testing.parametrize("dtype", [np.float64, np.float32, np.uint8, bool])
|
|
@testing.parametrize("preserve_range", [False, True])
|
|
def test_radon_center(shape, circle, dtype, preserve_range):
|
|
check_radon_center(shape, circle, dtype, preserve_range)
|
|
|
|
|
|
@testing.parametrize("shape", [(32, 16), (33, 17)])
|
|
@testing.parametrize("circle", [False])
|
|
@testing.parametrize("dtype", [np.float64, np.float32, np.uint8, bool])
|
|
@testing.parametrize("preserve_range", [False, True])
|
|
def test_radon_center_rectangular(shape, circle, dtype, preserve_range):
|
|
check_radon_center(shape, circle, dtype, preserve_range)
|
|
|
|
|
|
def check_iradon_center(size, theta, circle):
|
|
debug = False
|
|
# 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=float)
|
|
sinogram[size // 2, 0] = 1.
|
|
else:
|
|
diagonal = int(np.ceil(np.sqrt(2) * size))
|
|
sinogram = np.zeros((diagonal, 1), dtype=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()
|
|
|
|
assert np.allclose(reconstruction, reconstruction_opposite)
|
|
|
|
|
|
sizes_for_test_iradon_center = [16, 17]
|
|
thetas_for_test_iradon_center = [0, 90]
|
|
circles_for_test_iradon_center = [False, True]
|
|
|
|
|
|
@testing.parametrize("size, theta, circle",
|
|
itertools.product(sizes_for_test_iradon_center,
|
|
thetas_for_test_iradon_center,
|
|
circles_for_test_iradon_center))
|
|
def test_iradon_center(size, theta, circle):
|
|
check_iradon_center(size, theta, circle)
|
|
|
|
|
|
def check_radon_iradon(interpolation_type, filter_type):
|
|
debug = False
|
|
image = PHANTOM
|
|
reconstructed = iradon(radon(image, circle=False), filter_name=filter_type,
|
|
interpolation=interpolation_type, circle=False)
|
|
delta = np.mean(np.abs(image - reconstructed))
|
|
print('\n\tmean error:', delta)
|
|
if debug:
|
|
_debug_plot(image, reconstructed)
|
|
if filter_type in ('ramp', 'shepp-logan'):
|
|
if interpolation_type == 'nearest':
|
|
allowed_delta = 0.03
|
|
else:
|
|
allowed_delta = 0.025
|
|
else:
|
|
allowed_delta = 0.05
|
|
assert delta < allowed_delta
|
|
|
|
|
|
filter_types = ["ramp", "shepp-logan", "cosine", "hamming", "hann"]
|
|
interpolation_types = ['linear', 'nearest']
|
|
radon_iradon_inputs = list(itertools.product(interpolation_types,
|
|
filter_types))
|
|
# cubic interpolation is slow; only run one test for it
|
|
radon_iradon_inputs.append(('cubic', 'shepp-logan'))
|
|
|
|
|
|
@testing.parametrize("interpolation_type, filter_type",
|
|
radon_iradon_inputs)
|
|
def test_radon_iradon(interpolation_type, filter_type):
|
|
check_radon_iradon(interpolation_type, filter_type)
|
|
|
|
|
|
@pytest.mark.parametrize("filter_type", filter_types)
|
|
def test_iradon_new_signature(filter_type):
|
|
image = PHANTOM
|
|
sinogram = radon(image, circle=False)
|
|
with pytest.warns(FutureWarning):
|
|
assert np.array_equal(iradon(sinogram, filter=filter_type),
|
|
iradon(sinogram, filter_name=filter_type))
|
|
|
|
|
|
def test_iradon_angles():
|
|
"""
|
|
Test with different number of projections
|
|
"""
|
|
size = 100
|
|
# Synthetic data
|
|
image = np.tri(size) + np.tri(size)[::-1]
|
|
# Large number of projections: a good quality is expected
|
|
nb_angles = 200
|
|
theta = np.linspace(0, 180, nb_angles, endpoint=False)
|
|
radon_image_200 = radon(image, theta=theta, circle=False)
|
|
reconstructed = iradon(radon_image_200, circle=False)
|
|
delta_200 = np.mean(abs(_rescale_intensity(image) -
|
|
_rescale_intensity(reconstructed)))
|
|
assert delta_200 < 0.03
|
|
# Lower number of projections
|
|
nb_angles = 80
|
|
radon_image_80 = radon(image, theta=theta, circle=False)
|
|
# Test whether the sum of all projections is approximately the same
|
|
s = radon_image_80.sum(axis=0)
|
|
assert np.allclose(s, s[0], rtol=0.01)
|
|
reconstructed = iradon(radon_image_80, circle=False)
|
|
delta_80 = np.mean(abs(image / np.max(image) -
|
|
reconstructed / np.max(reconstructed)))
|
|
# Loss of quality when the number of projections is reduced
|
|
assert delta_80 > delta_200
|
|
|
|
|
|
def check_radon_iradon_minimal(shape, slices):
|
|
debug = False
|
|
theta = np.arange(180)
|
|
image = np.zeros(shape, dtype=float)
|
|
image[slices] = 1.
|
|
sinogram = radon(image, theta, circle=False)
|
|
reconstructed = iradon(sinogram, theta, circle=False)
|
|
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))
|
|
|
|
|
|
shapes = [(3, 3), (4, 4), (5, 5)]
|
|
|
|
|
|
def generate_test_data_for_radon_iradon_minimal(shapes):
|
|
def shape2coordinates(shape):
|
|
c0, c1 = shape[0] // 2, shape[1] // 2
|
|
coordinates = itertools.product((c0 - 1, c0, c0 + 1),
|
|
(c1 - 1, c1, c1 + 1))
|
|
return coordinates
|
|
|
|
def shape2shapeandcoordinates(shape):
|
|
return itertools.product([shape], shape2coordinates(shape))
|
|
|
|
return itertools.chain.from_iterable([shape2shapeandcoordinates(shape)
|
|
for shape in shapes])
|
|
|
|
|
|
@testing.parametrize("shape, coordinate",
|
|
generate_test_data_for_radon_iradon_minimal(shapes))
|
|
def test_radon_iradon_minimal(shape, coordinate):
|
|
check_radon_iradon_minimal(shape, coordinate)
|
|
|
|
|
|
def test_reconstruct_with_wrong_angles():
|
|
a = np.zeros((3, 3))
|
|
p = radon(a, theta=[0, 1, 2], circle=False)
|
|
iradon(p, theta=[0, 1, 2], circle=False)
|
|
with testing.raises(ValueError):
|
|
iradon(p, theta=[0, 1, 2, 3])
|
|
|
|
|
|
def _random_circle(shape):
|
|
# Synthetic random data, zero outside reconstruction circle
|
|
np.random.seed(98312871)
|
|
image = np.random.rand(*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.
|
|
return image
|
|
|
|
|
|
def test_radon_circle():
|
|
a = np.ones((10, 10))
|
|
with expected_warnings(['reconstruction circle']):
|
|
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_intensity(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
|
|
image = _random_circle(shape)
|
|
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 < 3.2e-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)
|
|
|
|
def argmax_shape(a):
|
|
return 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))
|
|
|
|
|
|
@testing.parametrize("size", (50, 51))
|
|
def test_sinogram_circle_to_square(size):
|
|
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)
|
|
|
|
|
|
# if adding more shapes to test data, you might want to look at commit d0f2bac3f
|
|
shapes_radon_iradon_circle = ((61, 79), )
|
|
interpolations = ('nearest', 'linear')
|
|
output_sizes = (None,
|
|
min(shapes_radon_iradon_circle[0]),
|
|
max(shapes_radon_iradon_circle[0]),
|
|
97)
|
|
|
|
|
|
@testing.parametrize("shape, interpolation, output_size",
|
|
itertools.product(shapes_radon_iradon_circle,
|
|
interpolations, output_sizes))
|
|
def test_radon_iradon_circle(shape, interpolation, output_size):
|
|
check_radon_iradon_circle(interpolation, shape, output_size)
|
|
|
|
|
|
def test_order_angles_golden_ratio():
|
|
from skimage.transform.radon_transform import order_angles_golden_ratio
|
|
np.random.seed(1231)
|
|
lengths = [1, 4, 10, 180]
|
|
for l in lengths:
|
|
theta_ordered = np.linspace(0, 180, l, endpoint=False)
|
|
theta_random = np.random.uniform(0, 180, l)
|
|
for theta in (theta_random, theta_ordered):
|
|
indices = [x for x in order_angles_golden_ratio(theta)]
|
|
# no duplicate indices allowed
|
|
assert len(indices) == len(set(indices))
|
|
|
|
|
|
@test_parallel()
|
|
def test_iradon_sart():
|
|
debug = False
|
|
|
|
image = rescale(PHANTOM, 0.8, mode='reflect',
|
|
multichannel=False, anti_aliasing=False)
|
|
theta_ordered = np.linspace(0., 180., image.shape[0], endpoint=False)
|
|
theta_missing_wedge = np.linspace(0., 150., image.shape[0], endpoint=True)
|
|
for theta, error_factor in ((theta_ordered, 1.),
|
|
(theta_missing_wedge, 2.)):
|
|
sinogram = radon(image, theta, circle=True)
|
|
reconstructed = iradon_sart(sinogram, theta)
|
|
|
|
if debug:
|
|
from matplotlib import pyplot as plt
|
|
plt.figure()
|
|
plt.subplot(221)
|
|
plt.imshow(image, interpolation='nearest')
|
|
plt.subplot(222)
|
|
plt.imshow(sinogram, interpolation='nearest')
|
|
plt.subplot(223)
|
|
plt.imshow(reconstructed, interpolation='nearest')
|
|
plt.subplot(224)
|
|
plt.imshow(reconstructed - image, interpolation='nearest')
|
|
plt.show()
|
|
|
|
delta = np.mean(np.abs(reconstructed - image))
|
|
print('delta (1 iteration) =', delta)
|
|
assert delta < 0.02 * error_factor
|
|
reconstructed = iradon_sart(sinogram, theta, reconstructed)
|
|
delta = np.mean(np.abs(reconstructed - image))
|
|
print('delta (2 iterations) =', delta)
|
|
assert delta < 0.014 * error_factor
|
|
reconstructed = iradon_sart(sinogram, theta, clip=(0, 1))
|
|
delta = np.mean(np.abs(reconstructed - image))
|
|
print('delta (1 iteration, clip) =', delta)
|
|
assert delta < 0.018 * error_factor
|
|
|
|
np.random.seed(1239867)
|
|
shifts = np.random.uniform(-3, 3, sinogram.shape[1])
|
|
x = np.arange(sinogram.shape[0])
|
|
sinogram_shifted = np.vstack([np.interp(x + shifts[i], x,
|
|
sinogram[:, i])
|
|
for i in range(sinogram.shape[1])]).T
|
|
reconstructed = iradon_sart(sinogram_shifted, theta,
|
|
projection_shifts=shifts)
|
|
if debug:
|
|
from matplotlib import pyplot as plt
|
|
plt.figure()
|
|
plt.subplot(221)
|
|
plt.imshow(image, interpolation='nearest')
|
|
plt.subplot(222)
|
|
plt.imshow(sinogram_shifted, interpolation='nearest')
|
|
plt.subplot(223)
|
|
plt.imshow(reconstructed, interpolation='nearest')
|
|
plt.subplot(224)
|
|
plt.imshow(reconstructed - image, interpolation='nearest')
|
|
plt.show()
|
|
|
|
delta = np.mean(np.abs(reconstructed - image))
|
|
print('delta (1 iteration, shifted sinogram) =', delta)
|
|
assert delta < 0.022 * error_factor
|
|
|
|
|
|
@pytest.mark.parametrize("preserve_range", [True, False])
|
|
def test_iradon_dtype(preserve_range):
|
|
sinogram = np.zeros((16, 1), dtype=int)
|
|
sinogram[8, 0] = 1.
|
|
sinogram64 = sinogram.astype('float64')
|
|
sinogram32 = sinogram.astype('float32')
|
|
|
|
assert iradon(sinogram, theta=[0],
|
|
preserve_range=preserve_range).dtype == 'float64'
|
|
assert iradon(sinogram64, theta=[0],
|
|
preserve_range=preserve_range).dtype == sinogram64.dtype
|
|
assert iradon(sinogram32, theta=[0],
|
|
preserve_range=preserve_range).dtype == sinogram32.dtype
|
|
|
|
|
|
def test_radon_dtype():
|
|
img = convert_to_float(PHANTOM, False)
|
|
img32 = img.astype(np.float32)
|
|
|
|
assert radon(img).dtype == img.dtype
|
|
assert radon(img32).dtype == img32.dtype
|
|
|
|
|
|
@pytest.mark.parametrize("dtype", [np.float32, np.float64])
|
|
def test_iradon_sart_dtype(dtype):
|
|
sinogram = np.zeros((16, 1), dtype=int)
|
|
sinogram[8, 0] = 1.
|
|
sinogram64 = sinogram.astype('float64')
|
|
sinogram32 = sinogram.astype('float32')
|
|
|
|
with expected_warnings(['Input data is cast to float']):
|
|
assert iradon_sart(sinogram, theta=[0]).dtype == 'float64'
|
|
|
|
assert iradon_sart(sinogram64, theta=[0]).dtype == sinogram64.dtype
|
|
assert iradon_sart(sinogram32, theta=[0]).dtype == sinogram32.dtype
|
|
|
|
assert iradon_sart(sinogram, theta=[0], dtype=dtype).dtype == dtype
|
|
assert iradon_sart(sinogram32, theta=[0], dtype=dtype).dtype == dtype
|
|
assert iradon_sart(sinogram64, theta=[0], dtype=dtype).dtype == dtype
|
|
|
|
|
|
def test_iradon_sart_wrong_dtype():
|
|
sinogram = np.zeros((16, 1))
|
|
|
|
with testing.raises(ValueError):
|
|
iradon_sart(sinogram, dtype=int)
|