wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-01 07:27:48 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #596 from josteinbf/radon-projection-center

Browse Source 
Radon transform: Improve tests and address bugs in definition of image and projection center.
This commit is contained in:
Stefan van der Walt
2013-07-04 07:24:05 -07:00
parent 647aacb316 6d2f082c11
commit bcc249263f
2 changed files with 233 additions and 124 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/transform/radon_transform.py
+39 -46
View File
@@ -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)
skimage/transform/tests/test_radon_transform.py
+194 -78
View File
@@ -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__":