wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-06-29 14:56:02 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #654 from josteinbf/radon-square-center

Browse Source 
Define rotation axis center for Radon transforms with inverses.
This commit is contained in:
Stefan van der Walt
2013-08-06 02:04:32 -07:00
parent 697d7dddf9 d5b72f91ab
commit be14f0b28c
2 changed files with 69 additions and 52 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/transform/radon_transform.py
+65 -49
View File
@@ -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
@@ -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
------
@@ -61,10 +64,11 @@ 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')
# Crop image to make it square
slices = []
for d in (0, 1):
if image.shape[d] > min(image.shape):
@@ -76,24 +80,25 @@ 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]
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
# 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):
@@ -105,14 +110,17 @@ 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):
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)
@@ -128,7 +136,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)).
@@ -148,8 +159,10 @@ def iradon(radon_image, theta=None, output_size=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
-----
@@ -190,39 +203,35 @@ 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
fourier_filter[:] = 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
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
[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
@@ -239,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))
@@ -258,9 +267,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
-----
@@ -314,7 +325,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)).
@@ -336,8 +350,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
-----
@@ -358,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)
skimage/transform/tests/test_radon_transform.py
+4 -3
View File
@@ -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
@@ -210,7 +211,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 +237,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 +285,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)