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

Merge pull request #567 from josteinbf/radon-reconstruction-circle

Browse Source 
ENH: Option to use a reconstruction circle in transform.radon/iradon.
This commit is contained in:
Stefan van der Walt
2013-05-27 13:14:10 -07:00
parent ba354dde4d d5629b1265
commit 7ced6625d8
3 changed files with 147 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CONTRIBUTORS.txt
+3
View File
@@ -139,3 +139,6 @@
- Xavier Moles Lopez
Color separation (color deconvolution) for several stainings.
- Jostein Bø Fløystad
Reconstruction circle mode for Radon transform
skimage/transform/radon_transform.py
+75 -20
View File
@@ -20,7 +20,7 @@ from ._warps_cy import _warp_fast
__all__ = ["radon", "iradon"]
def radon(image, theta=None):
def radon(image, theta=None, circle=False):
"""
Calculates the radon transform of an image given specified
projection angles.
@@ -31,31 +31,61 @@ def radon(image, theta=None):
Input image.
theta : array_like, dtype=float, optional (default np.arange(180))
Projection angles (in degrees).
circle : boolean, optional
Assume image is zero outside the inscribed circle, making the
width of each projection (the first dimension of the sinogram)
equal to ``min(image.shape)``.
Returns
-------
output : ndarray
Radon transform (sinogram).
Raises
------
ValueError
If called with ``circle=True`` and ``image != 0`` outside the inscribed
circle
"""
if image.ndim != 2:
raise ValueError('The input image must be 2-D')
if theta is None:
theta = np.arange(180)
height, width = image.shape
diagonal = np.sqrt(height**2 + width**2)
heightpad = np.ceil(diagonal - height)
widthpad = np.ceil(diagonal - width)
padded_image = np.zeros((int(height + heightpad),
int(width + widthpad)))
y0, y1 = int(np.ceil(heightpad / 2)), \
int((np.ceil(heightpad / 2) + height))
x0, x1 = int((np.ceil(widthpad / 2))), \
int((np.ceil(widthpad / 2) + width))
if circle:
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
if not np.all(reconstruction_circle | (image == 0)):
raise ValueError('Image must be zero outside the reconstruction'
' circle')
slices = []
for d in (0, 1):
if image.shape[d] > min(image.shape):
excess = image.shape[d] - min(image.shape)
slices.append(slice(int(np.ceil(excess / 2)),
int(np.ceil(excess / 2)
+ min(image.shape))))
else:
slices.append(slice(None))
slices = tuple(slices)
padded_image = image[slices]
out = np.zeros((min(padded_image.shape), len(theta)))
else:
height, width = image.shape
diagonal = np.sqrt(height**2 + width**2)
heightpad = np.ceil(diagonal - height)
widthpad = 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))
padded_image[y0:y1, x0:x1] = image
out = np.zeros((max(padded_image.shape), len(theta)))
padded_image[y0:y1, x0:x1] = image
out = np.zeros((max(padded_image.shape), len(theta)))
h, w = padded_image.shape
dh, dw = h // 2, w // 2
@@ -86,7 +116,7 @@ def radon(image, theta=None):
def iradon(radon_image, theta=None, output_size=None,
filter="ramp", interpolation="linear"):
filter="ramp", interpolation="linear", circle=False):
"""
Inverse radon transform.
@@ -110,6 +140,10 @@ def iradon(radon_image, theta=None, output_size=None,
interpolation : str, optional (default linear)
Interpolation method used in reconstruction.
Methods available: nearest, linear.
circle : boolean, optional
Assume the reconstructed image is zero outside the inscribed circle.
Also changes the default output_size to match the behaviour of
``radon`` called with ``circle=True``.
Returns
-------
@@ -139,7 +173,19 @@ def iradon(radon_image, theta=None, output_size=None,
th = (np.pi / 180.0) * theta
# if output size not specified, estimate from input radon image
if not output_size:
output_size = int(np.floor(np.sqrt((radon_image.shape[0])**2 / 2.0)))
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
n = radon_image.shape[0]
img = radon_image.copy()
@@ -163,7 +209,7 @@ def iradon(radon_image, theta=None, output_size=None,
f[1:] = f[1:] * (0.54 + 0.46 * np.cos(w[1:]))
elif filter == "hann":
f[1:] = f[1:] * (1 + np.cos(w[1:])) / 2
elif filter == None:
elif filter is None:
f[1:] = 1
else:
raise ValueError("Unknown filter: %s" % filter)
@@ -185,12 +231,19 @@ def iradon(radon_image, theta=None, output_size=None,
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
if interpolation == "nearest":
for i in range(len(theta)):
k = np.round(mid_index + xpr * np.sin(th[i]) - ypr * np.cos(th[i]))
reconstructed += radon_filtered[
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)):
@@ -199,9 +252,11 @@ def iradon(radon_image, theta=None, output_size=None,
b = mid_index + a
b0 = ((((b + 1 > 0) & (b + 1 < n)) * (b + 1)) - 1).astype(np.int)
b1 = ((((b > 0) & (b < n)) * b) - 1).astype(np.int)
reconstructed += (t - a) * radon_filtered[b0, i] + \
(a - t + 1) * radon_filtered[b1, i]
backprojected = (t - a) * radon_filtered[b0, i] + \
(a - t + 1) * radon_filtered[b1, i]
if circle:
backprojected[~reconstruction_circle] = 0.
reconstructed += backprojected
else:
raise ValueError("Unknown interpolation: %s" % interpolation)
skimage/transform/tests/test_radon_transform.py
+69
View File
@@ -1,7 +1,9 @@
from __future__ import print_function
from __future__ import division
import numpy as np
from numpy.testing import *
import itertools
from skimage.transform import *
@@ -104,5 +106,72 @@ def test_reconstruct_with_wrong_angles():
assert_raises(ValueError, iradon, p, theta=[0, 1, 2, 3])
def test_radon_circle():
a = np.ones((10, 10))
assert_raises(ValueError, 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(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
np.random.seed(98312871)
image = np.random.rand(*shape)
image[r >= radius] = 0.
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 < 3e-3)
def test_radon_iradon_circle():
shape = (61, 79)
# Synthetic random data, zero outside reconstruction circle
image = np.random.rand(*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)
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.
# 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)
if __name__ == "__main__":
run_module_suite()