mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-29 11:26:57 +08:00
Merge branch 'stefan-radon'
This commit is contained in:
Emmanuelle Gouillart
@@ -0,0 +1,116 @@
|
||||
***************
|
||||
Radon transform
|
||||
***************
|
||||
|
||||
The radon transform is a technique widely used in tomography, where you
|
||||
reconstruct an object from its different projections. A projection for example
|
||||
the scattering data obtained as the output of a tomographic scan.
|
||||
|
||||
For more information:
|
||||
http://en.wikipedia.org/wiki/Radon_transform
|
||||
http://www.clear.rice.edu/elec431/projects96/DSP/bpanalysis.html
|
||||
|
||||
Forward transform
|
||||
=================
|
||||
|
||||
First we load the Schepp-Logan phantom, a classic test image representing a
|
||||
tomographic scan.
|
||||
|
||||
.. ipython::
|
||||
|
||||
In [1]: from scikits.image.io import imread
|
||||
|
||||
In [1]: from scikits.image import data_dir
|
||||
|
||||
In [2]: from scikits.image.transform import radon, iradon
|
||||
|
||||
In [3]: from scikits.image.color import rgb2gray
|
||||
|
||||
In [4]: import matplotlib.pyplot as plt
|
||||
|
||||
In [5]: import matplotlib.cm as cm
|
||||
|
||||
In [6]: image = rgb2gray(imread(data_dir + "/phantom.png"))
|
||||
|
||||
In [7]: plt.title("original image");
|
||||
|
||||
In [8]: plt.imshow(image, cmap=cm.Greys_r)
|
||||
|
||||
@savefig radon_original_image.png width=4in
|
||||
In [9]: plt.show()
|
||||
|
||||
|
||||
Let us illustrate the transform by looking at projections taken at specific
|
||||
angles.
|
||||
|
||||
.. ipython::
|
||||
|
||||
In [10]: projections = radon(image, theta=[0, 45, 90])
|
||||
|
||||
In [11]: plt.plot(projections);
|
||||
|
||||
In [12]: plt.title("radon projections");
|
||||
|
||||
In [13]: plt.xlabel("projection axis");
|
||||
|
||||
In [14]: plt.ylabel("intensity");
|
||||
|
||||
@savefig radon_projection_plot1.png width=4in
|
||||
In [15]: plt.show()
|
||||
|
||||
We are going to reconstruct an image from 180 of these projections (the
|
||||
default).
|
||||
|
||||
.. ipython::
|
||||
|
||||
In [16]: projections = radon(image)
|
||||
|
||||
In [17]: plt.figure()
|
||||
|
||||
In [18]: plt.title("radon projections");
|
||||
|
||||
In [19]: plt.xlabel("projection axis");
|
||||
|
||||
In [20]: plt.ylabel("intensity");
|
||||
|
||||
In [21]: plt.plot(projections)
|
||||
|
||||
@savefig radon_projection_plot2.png width=4in
|
||||
In [22]: plt.show()
|
||||
|
||||
|
||||
We have now constructed various projections, line integrals of an image, at
|
||||
specific angles. This image is called a sinogram.
|
||||
|
||||
.. ipython::
|
||||
|
||||
In [23]: plt.figure()
|
||||
|
||||
In [24]: plt.title("sinogram");
|
||||
|
||||
In [25]: plt.xlabel("projection axis");
|
||||
|
||||
In [26]: plt.ylabel("intensity");
|
||||
|
||||
In [27]: plt.imshow(projections)
|
||||
|
||||
@savefig radon_sinogram.png width=4in
|
||||
In [28]: plt.show()
|
||||
|
||||
|
||||
Inverse transform
|
||||
=================
|
||||
To reconstruct the image from this sinogram, we apply the inverse transform.
|
||||
|
||||
.. ipython::
|
||||
|
||||
In [29]: reconstruction = iradon(projections)
|
||||
|
||||
In [30]: plt.title("reconstructed image");
|
||||
|
||||
In [31]: plt.imshow(reconstruction, cmap=cm.Greys_r)
|
||||
|
||||
@savefig radon_reconstructed_image.png width=4in
|
||||
In [32]: plt.show()
|
||||
|
||||
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 3.3 KiB |
@@ -1,4 +1,6 @@
|
||||
from .hough_transform import *
|
||||
from .radon_transform import *
|
||||
from .finite_radon_transform import *
|
||||
from .project import *
|
||||
from ._project import homography as fast_homography
|
||||
from .integral import *
|
||||
|
||||
@@ -0,0 +1,209 @@
|
||||
#cython: cdivison=True boundscheck=False
|
||||
|
||||
__all__ = ['homography']
|
||||
|
||||
cimport cython
|
||||
cimport numpy as np
|
||||
|
||||
import numpy as np
|
||||
import cython
|
||||
|
||||
from cython.operator import dereference
|
||||
|
||||
np.import_array()
|
||||
|
||||
cdef extern from "math.h":
|
||||
double floor(double)
|
||||
double fmod(double, double)
|
||||
|
||||
cdef double get_pixel(double *image, int rows, int cols,
|
||||
int r, int c, char mode, double cval=0):
|
||||
"""Get a pixel from the image, taking wrapping mode into consideration.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : *double
|
||||
Input image.
|
||||
rows, cols : int
|
||||
Dimensions of image.
|
||||
r, c : int
|
||||
Position at which to get the pixel.
|
||||
mode : {'C', 'W', 'M'}
|
||||
Wrapping mode. Constant, Wrap or Mirror.
|
||||
cval : double
|
||||
Constant value to use for mode constant.
|
||||
|
||||
"""
|
||||
if mode == 'C':
|
||||
if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1):
|
||||
return cval
|
||||
else:
|
||||
return image[r * cols + c]
|
||||
else:
|
||||
return image[coord_map(rows, r, mode) * cols +
|
||||
coord_map(cols, c, mode)]
|
||||
|
||||
cdef int coord_map(int dim, int coord, char mode):
|
||||
"""
|
||||
Wrap a coordinate, according to a given dimension and mode.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
dim : int
|
||||
Maximum coordinate.
|
||||
coord : int
|
||||
Coord provided by user. May be < 0 or > dim.
|
||||
mode : {'W', 'M'}
|
||||
Whether to wrap or mirror the coordinate if it
|
||||
falls outside [0, dim).
|
||||
|
||||
"""
|
||||
dim = dim - 1
|
||||
if mode == 'M': # mirror
|
||||
if (coord < 0):
|
||||
# How many times times does the coordinate wrap?
|
||||
if (<int>(-coord / dim) % 2 != 0):
|
||||
return dim - <int>(-coord % dim)
|
||||
else:
|
||||
return <int>(-coord % dim)
|
||||
elif (coord > dim):
|
||||
if (<int>(coord / dim) % 2 != 0):
|
||||
return <int>(dim - (coord % dim))
|
||||
else:
|
||||
return <int>(coord % dim)
|
||||
elif mode == 'W': # wrap
|
||||
if (coord < 0):
|
||||
return <int>(dim - (-coord % dim))
|
||||
elif (coord > dim):
|
||||
return <int>(coord % dim)
|
||||
|
||||
return coord
|
||||
|
||||
cdef tf(double x, double y, double* H, double *x_, double *y_):
|
||||
"""Apply a homography to a coordinate.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x, y : double
|
||||
Input coordinate.
|
||||
H : (3,3) *double
|
||||
Transformation matrix.
|
||||
x_, y_ : *double
|
||||
Output coordinate.
|
||||
|
||||
"""
|
||||
cdef double xx, yy, zz
|
||||
|
||||
xx = H[0] * x + H[1] * y + H[2]
|
||||
yy = H[3] * x + H[4] * y + H[5]
|
||||
zz = H[6] * x + H[7] * y + H[8]
|
||||
|
||||
xx = xx / zz
|
||||
yy = yy / zz
|
||||
|
||||
x_[0] = xx
|
||||
y_[0] = yy
|
||||
|
||||
@cython.boundscheck(False)
|
||||
def homography(np.ndarray image, np.ndarray H, output_shape=None,
|
||||
mode='constant', double cval=0):
|
||||
"""
|
||||
Projective transformation (homography).
|
||||
|
||||
Perform a projective transformation (homography) of a
|
||||
floating point image, using bi-linear interpolation.
|
||||
|
||||
For each pixel, given its homogeneous coordinate :math:`\mathbf{x}
|
||||
= [x, y, 1]^T`, its target position is calculated by multiplying
|
||||
with the given matrix, :math:`H`, to give :math:`H \mathbf{x}`.
|
||||
E.g., to rotate by theta degrees clockwise, the matrix should be
|
||||
|
||||
::
|
||||
|
||||
[[cos(theta) -sin(theta) 0]
|
||||
[sin(theta) cos(theta) 0]
|
||||
[0 0 1]]
|
||||
|
||||
or, to translate x by 10 and y by 20,
|
||||
|
||||
::
|
||||
|
||||
[[1 0 10]
|
||||
[0 1 20]
|
||||
[0 0 1 ]].
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : 2-D array
|
||||
Input image.
|
||||
H : array of shape ``(3, 3)``
|
||||
Transformation matrix H that defines the homography.
|
||||
output_shape : tuple (rows, cols)
|
||||
Shape of the output image generated.
|
||||
order : int
|
||||
Order of splines used in interpolation.
|
||||
mode : {'constant', 'mirror', 'wrap'}
|
||||
How to handle values outside the image borders.
|
||||
cval : string
|
||||
Used in conjunction with mode 'C' (constant), the value
|
||||
outside the image boundaries.
|
||||
|
||||
"""
|
||||
|
||||
cdef np.ndarray[dtype=np.double_t, ndim=2] img = \
|
||||
np.asarray(image, dtype=np.double)
|
||||
cdef np.ndarray[dtype=np.double_t, ndim=2] M = \
|
||||
np.ascontiguousarray(np.linalg.inv(H))
|
||||
|
||||
if mode not in ('constant', 'wrap', 'mirror'):
|
||||
raise ValueError("Invalid mode specified. Please use "
|
||||
"`constant`, `wrap` or `mirror`.")
|
||||
if mode == 'constant':
|
||||
mode_c = ord('C')
|
||||
elif mode == 'wrap':
|
||||
mode_c = ord('W')
|
||||
elif mode == 'mirror':
|
||||
mode_c = ord('M')
|
||||
|
||||
cdef int out_r, out_c, columns, rows
|
||||
if output_shape is None:
|
||||
out_r = img.shape[0]
|
||||
out_c = img.shape[1]
|
||||
else:
|
||||
out_r = output_shape[0]
|
||||
out_c = output_shape[1]
|
||||
|
||||
rows = img.shape[0]
|
||||
columns = img.shape[1]
|
||||
|
||||
cdef np.ndarray[dtype=np.double_t, ndim=2] out = \
|
||||
np.zeros((out_r, out_c), dtype=np.double)
|
||||
|
||||
cdef int tfr, tfc, r_int, c_int
|
||||
cdef double y0, y1, y2, y3
|
||||
cdef double r, c, z, t, u
|
||||
|
||||
for tfr in range(out_r):
|
||||
for tfc in range(out_c):
|
||||
tf(tfc, tfr, <double*>M.data, &c, &r)
|
||||
r_int = <int>floor(r)
|
||||
c_int = <int>floor(c)
|
||||
|
||||
t = r - r_int
|
||||
u = c - c_int
|
||||
|
||||
y0 = get_pixel(<double*>img.data, rows, columns,
|
||||
r_int, c_int, mode_c)
|
||||
y1 = get_pixel(<double*>img.data, rows, columns,
|
||||
r_int + 1, c_int, mode_c)
|
||||
y2 = get_pixel(<double*>img.data, rows, columns,
|
||||
r_int + 1, c_int + 1, mode_c)
|
||||
y3 = get_pixel(<double*>img.data, rows, columns,
|
||||
r_int, c_int + 1, mode_c)
|
||||
|
||||
out[tfr, tfc] = \
|
||||
(1 - t) * (1 - u) * y0 + \
|
||||
t * (1 - u) * y1 + \
|
||||
t * u * y2 + (1 - t) * u * y3;
|
||||
|
||||
return out
|
||||
@@ -0,0 +1,191 @@
|
||||
"""
|
||||
radon.py - Radon and inverse radon transforms
|
||||
|
||||
Based on code of Justin K. Romberg
|
||||
(http://www.clear.rice.edu/elec431/projects96/DSP/bpanalysis.html)
|
||||
J. Gillam and Chris Griffin.
|
||||
|
||||
References:
|
||||
-B.R. Ramesh, N. Srinivasa, K. Rajgopal, "An Algorithm for Computing
|
||||
the Discrete Radon Transform With Some Applications", Proceedings of
|
||||
the Fourth IEEE Region 10 International Conference, TENCON '89, 1989.
|
||||
-A. C. Kak, Malcolm Slaney, "Principles of Computerized Tomographic
|
||||
Imaging", IEEE Press 1988.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
from scipy.fftpack import fftshift, fft, ifft
|
||||
from ._project import homography
|
||||
|
||||
def radon(image, theta=None):
|
||||
"""
|
||||
Calculates the radon transform of an image given specified
|
||||
projection angles.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
image : array_like, dtype=float
|
||||
Input image.
|
||||
theta : array_like, dtype=float, optional (default np.arange(180))
|
||||
Projection angles (in degrees).
|
||||
|
||||
Returns
|
||||
-------
|
||||
output : ndarray
|
||||
Radon transform (sinogram).
|
||||
|
||||
"""
|
||||
if image.ndim != 2:
|
||||
raise ValueError('The input image must be 2-D')
|
||||
if theta == None:
|
||||
theta = np.arange(180)
|
||||
height, width = image.shape
|
||||
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[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 = 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"):
|
||||
"""
|
||||
Inverse radon transform.
|
||||
|
||||
Reconstruct an image from the radon transform, using the filtered
|
||||
back projection algorithm.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
radon_image : array_like, dtype=float
|
||||
Image containing radon transform (sinogram). Each column of
|
||||
the image corresponds to a projection along a different angle.
|
||||
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 nxm)
|
||||
output_size : int
|
||||
Number of rows and columns in the reconstruction.
|
||||
filter : str, optional (default ramp)
|
||||
Filter used in frequency domain filtering. Ramp filter used by default.
|
||||
Filters available: ramp, shepp-logan, cosine, hamming, hann
|
||||
Assign None to use no filter.
|
||||
interpolation : str, optional (default linear)
|
||||
Interpolation method used in reconstruction.
|
||||
Methods available: nearest, linear.
|
||||
|
||||
Returns
|
||||
-------
|
||||
output : ndarray
|
||||
Reconstructed image.
|
||||
|
||||
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. This algorithm is called filtered back projection.
|
||||
|
||||
"""
|
||||
if radon_image.ndim != 2:
|
||||
raise ValueError('The input image must be 2-D')
|
||||
if theta == None:
|
||||
m, n = radon_image.shape
|
||||
theta = np.linspace(0, 180, n, endpoint=False)
|
||||
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)))
|
||||
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
|
||||
freqs = np.zeros((order, 1))
|
||||
|
||||
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
|
||||
if filter == "ramp":
|
||||
pass
|
||||
elif filter == "shepp-logan":
|
||||
f[1:] = f[1:] * np.sin(w[1:] / 2) / (w[1:] / 2)
|
||||
elif filter == "cosine":
|
||||
f[1:] = f[1:] * np.cos(w[1:] / 2)
|
||||
elif filter == "hamming":
|
||||
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:
|
||||
f[1:] = 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
|
||||
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))
|
||||
mid_index = np.ceil(n/2);
|
||||
x = output_size
|
||||
y = output_size
|
||||
[X, Y] = np.mgrid[0.0:x, 0.0:y]
|
||||
xpr = X - (output_size + 1.0) / 2.0
|
||||
ypr = Y - (output_size + 1.0) / 2.0
|
||||
|
||||
# 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]
|
||||
elif interpolation == "linear":
|
||||
for i in range(len(theta)):
|
||||
t = xpr*np.sin(th[i]) - ypr*np.cos(th[i])
|
||||
a = np.floor(t)
|
||||
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]
|
||||
else:
|
||||
raise ValueError("Unknown interpolation: %s" % interpolation)
|
||||
|
||||
return reconstructed * np.pi / (2 * len(th))
|
||||
@@ -15,10 +15,14 @@ def configuration(parent_package='', top_path=None):
|
||||
config.add_data_dir('tests')
|
||||
|
||||
cython(['_hough_transform.pyx'], working_path=base_path)
|
||||
cython(['_project.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('_hough_transform', sources=['_hough_transform.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
config.add_extension('_project', sources=['_project.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
return config
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
@@ -1,7 +1,9 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_array_almost_equal
|
||||
|
||||
from scikits.image.transform.project import _stackcopy, homography
|
||||
from scikits.image.transform.project import _stackcopy
|
||||
from scikits.image.transform import homography, fast_homography
|
||||
from scikits.image import data
|
||||
|
||||
def test_stackcopy():
|
||||
layers = 4
|
||||
@@ -19,3 +21,41 @@ def test_homography():
|
||||
[0, 0, 1]])
|
||||
x90 = homography(x, M, order=1)
|
||||
assert_array_almost_equal(x90, np.rot90(x))
|
||||
|
||||
def test_fast_homography():
|
||||
img = data.lena()
|
||||
img = img[:, :100]
|
||||
|
||||
theta = np.deg2rad(30)
|
||||
scale = 0.5
|
||||
tx, ty = 50, 50
|
||||
|
||||
H = np.eye(3)
|
||||
S = scale * np.sin(theta)
|
||||
C = scale * np.cos(theta)
|
||||
|
||||
H[:2, :2] = [[C, -S], [S, C]]
|
||||
H[:2, 2] = [tx, ty]
|
||||
|
||||
for mode in ('constant', 'mirror', 'wrap'):
|
||||
print 'Transform mode:', mode
|
||||
|
||||
p0 = homography(img, H, mode=mode, order=1)
|
||||
p1 = fast_homography(img, H, mode=mode)
|
||||
p1 = np.round(p1)
|
||||
|
||||
## import matplotlib.pyplot as plt
|
||||
## f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
|
||||
## ax0.imshow(img)
|
||||
## ax1.imshow(p0, cmap=plt.cm.gray)
|
||||
## ax2.imshow(p1, cmap=plt.cm.gray)
|
||||
## ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
|
||||
## plt.show()
|
||||
|
||||
d = np.mean(np.abs(p0 - p1))
|
||||
assert d < 0.2
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
from numpy.testing import run_module_suite
|
||||
run_module_suite()
|
||||
|
||||
@@ -0,0 +1,63 @@
|
||||
import numpy as np
|
||||
from numpy.testing import *
|
||||
from scikits.image.transform import *
|
||||
|
||||
def rescale(x):
|
||||
x = x.astype(float)
|
||||
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)
|
||||
|
||||
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.mean(abs(image - reconstructed))
|
||||
assert delta < 0.05
|
||||
|
||||
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
|
||||
radon_image_200 = radon(image, theta=np.linspace(0, 180, nb_angles,
|
||||
endpoint=False))
|
||||
reconstructed = iradon(radon_image_200)
|
||||
delta_200 = np.mean(abs(rescale(image) - rescale(reconstructed)))
|
||||
assert delta_200 < 0.03
|
||||
# Lower number of projections
|
||||
nb_angles = 80
|
||||
radon_image_80 = radon(image, theta=np.linspace(0, 180, nb_angles,
|
||||
endpoint=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)
|
||||
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
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
run_module_suite()
|
||||
|
||||
Reference in New Issue
Block a user