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

ENH: Use fast homography instead of PIL's rotate. Minor cleanup.

Browse Source
This commit is contained in:
Stefan van der Walt
2011-09-25 18:04:27 -07:00
parent ab19017962
commit b65af2727a
3 changed files with 76 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files
scikits/image/transform/__init__.py
+1 -1
View File
@@ -1,5 +1,5 @@
from .hough_transform import *
from .radon import *
from .radon_transform import *
from .finite_radon_transform import *
from .project import *
from ._project import homography as fast_homography
scikits/image/transform/radon_transform.py
+54 -24
View File
@@ -14,15 +14,13 @@ References:
"""
import numpy as np
from scipy.misc import imrotate
from scipy.interpolate import interp1d
from scipy.fftpack import fftshift, fft, ifft
import math
from ._project import homography
def radon(image, theta=None):
"""
Calculates the radon transform of an image given specified projection angles.
Calculates the radon transform of an image given specified
projection angles.
Parameters
----------
@@ -35,6 +33,7 @@ def radon(image, theta=None):
-------
output : ndarray
Radon transform.
"""
if image.ndim != 2:
raise ValueError('The input image must be 2-D')
@@ -44,25 +43,54 @@ def radon(image, theta=None):
diagonal = np.sqrt(height ** 2 + width ** 2)
heightpad = np.ceil(diagonal - height) + 2
widthpad = np.ceil(diagonal - width) + 2
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))
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))
padded_image[y0:y1, x0:x1] = image
out = np.zeros((max(padded_image.shape), len(theta)))
h, w = padded_image.shape
shift0 = np.array([[1, 0, -w/2.],
[0, 1, -h/2.],
[0, 0, 1]])
shift1 = np.array([[1, 0, w/2.],
[0, 1, h/2.],
[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],
[0, 0, 1]])
return shift1.dot(R).dot(shift0)
for i in range(len(theta)):
rotated = imrotate(padded_image, -theta[i])
rotated = homography(padded_image,
build_rotation(-theta[i]))
out[:,i] = rotated.sum(0)[::-1]
return out
def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolation="linear"):
def iradon(radon_image, theta=None, output_size=None,
filter="ramp", interpolation="linear"):
"""
Reconstructs an image from radon transformed data.
Inverse radon transform.
Reconstruct an image from the radon transform.
Parameters
----------
radon_image : array_like, dtype=float
Image containing radon transform.
Image containing radon transform (sinogram).
theta : array_like, dtype=float, optional (default np.arange(180))
Reconstruction angles (in degrees).
output_size : int
@@ -82,17 +110,19 @@ def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolati
Notes
-----
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.
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.
"""
if radon_image.ndim != 2:
raise ValueError('The input image must be 2-D')
if theta == None:
theta = np.arange(180)
th = (math.pi/180.0)*theta
th = (np.pi / 180.0) * theta
# if output size not specified, estimate from input radon image
if not output_size:
output_size = 2*np.floor(radon_image.shape[0] / (2 * np.sqrt(2)))
output_size = 2 * np.floor(radon_image.shape[0] / (2 * np.sqrt(2)))
n = radon_image.shape[0]
img = radon_image.copy()
@@ -101,11 +131,11 @@ def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolati
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
freqs = np.zeros((order, 1))
f = fftshift(abs(np.mgrid[-1:1:2 / order])).reshape(-1, 1)
w = 2 * math.pi * f
w = 2 * np.pi * f
# start from first element to avoid divide by zero
if filter == "ramp":
pass
@@ -139,8 +169,9 @@ def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolati
# 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[((((k > 0) & (k < n)) * k) - 1).astype(np.int), i]
k = np.round(mid_index + xpr * np.sin(th[i]) - ypr * np.cos(th[i]))
reconstructed += radon_filtered[
((((k > 0) & (k < n)) * k) - 1).astype(np.int), i]
elif interpolation == "linear":
for i in range(len(theta)):
t = xpr*np.sin(th[i]) - ypr*np.cos(th[i])
@@ -148,10 +179,9 @@ def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolati
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]
reconstructed += (t - a) * radon_filtered[b0, i] + \
(a - t + 1) * radon_filtered[b1, i]
else:
raise ValueError("Unknown interpolation: %s" % interpolation)
return reconstructed * math.pi / (2*len(th))
return reconstructed * np.pi / (2 * len(th))
scikits/image/transform/tests/test_radon_transform.py
+21 -4
View File
@@ -2,17 +2,34 @@ import numpy as np
from numpy.testing import *
from scikits.image.transform import *
def rescale(x):
x = x.astype(float, copy=True)
x -= x.min()
x /= x.max()
return x
def test_radon_iradon():
size = 100
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)
delta = np.sum(abs(image/np.max(image) - reconstructed/np.max(reconstructed)))/(size*size)
assert delta < 0.1
image = rescale(image)
reconstructed = rescale(reconstructed)
delta = np.mean(np.abs(image - reconstructed))
## 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
reconstructed = iradon(radon(image), filter="ramp", interpolation="nearest")
delta = np.sum(abs(image/np.max(image) - reconstructed/np.max(reconstructed)))/(size*size)
assert delta < 0.1
delta = np.mean(abs(image - reconstructed))
assert delta < 0.05
if __name__ == "__main__":