diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index e03271b8..ac11457c 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -144,3 +144,4 @@ - Jostein Bø Fløystad Reconstruction circle mode for Radon transform + Simultaneous Algebraic Reconstruction Technique for inverse Radon transform diff --git a/doc/examples/plot_radon_transform.py b/doc/examples/plot_radon_transform.py index b4b9721d..61175f8a 100644 --- a/doc/examples/plot_radon_transform.py +++ b/doc/examples/plot_radon_transform.py @@ -1,57 +1,200 @@ +# -*- coding: utf-8 -*- """ =============== Radon transform =============== -The radon transform is a technique widely used in tomography to -reconstruct an object from different projections. A projection is, for -example, the scattering data obtained as the output of a tomographic -scan. +In computed tomography, the tomography reconstruction problem is to obtain +a tomographic slice image from a set of projections [1]_. A projection is formed +by drawing a set of parallel rays through the 2D object of interest, assigning +the integral of the object's contrast along each ray to a single pixel in the +projection. A single projection of a 2D object is one dimensional. To +enable computed tomography reconstruction of the object, several projections +must be acquired, each of them corresponding to a different angle between the +rays with respect to the object. A collection of projections at several angles +is called a sinogram, which is a linear transform of the original image. -For more information see: +The inverse Radon transform is used in computed tomography to reconstruct +a 2D image from the measured projections (the sinogram). A practical, exact +implementation of the inverse Radon transform does not exist, but there are +several good approximate algorithms available. - - http://en.wikipedia.org/wiki/Radon_transform - - http://www.clear.rice.edu/elec431/projects96/DSP/bpanalysis.html +As the inverse Radon transform reconstructs the object from a set of +projections, the (forward) Radon transform can be used to simulate a +tomography experiment. -This script performs the radon transform, and reconstructs the -input image based on the resulting sinogram. +This script performs the Radon transform to simulate a tomography experiment +and reconstructs the input image based on the resulting sinogram formed by +the simulation. Two methods for performing the inverse Radon transform +and reconstructing the original image are compared: The Filtered Back +Projection (FBP) and the Simultaneous Algebraic Reconstruction +Technique (SART). +.. seealso:: + + - AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging", + http://www.slaney.org/pct/pct-toc.html + - http://en.wikipedia.org/wiki/Radon_transform + +The forward transform +===================== + +As our original image, we will use the Shepp-Logan phantom. When calculating +the Radon transform, we need to decide how many projection angles we wish +to use. As a rule of thumb, the number of projections should be about the +same as the number of pixels there are across the object (to see why this +is so, consider how many unknown pixel values must be determined in the +reconstruction process and compare this to the number of measurements +provided by the projections), and we follow that rule here. Below is the +original image and its Radon transform, often known as its _sinogram_: """ +from __future__ import print_function, division + +import numpy as np import matplotlib.pyplot as plt from skimage.io import imread from skimage import data_dir -from skimage.transform import radon, iradon, rescale - +from skimage.transform import radon, rescale image = imread(data_dir + "/phantom.png", as_grey=True) image = rescale(image, scale=0.4) -plt.figure(figsize=(8, 8.5)) +plt.figure(figsize=(8, 4.5)) -plt.subplot(221) +plt.subplot(121) plt.title("Original") plt.imshow(image, cmap=plt.cm.Greys_r) -plt.subplot(222) -projections = radon(image, theta=[0, 45, 90]) -plt.plot(projections) -plt.title("Projections at\n0, 45 and 90 degrees") -plt.xlabel("Projection axis") -plt.ylabel("Intensity") - -projections = radon(image) -plt.subplot(223) +theta = np.linspace(0., 180., max(image.shape), endpoint=True) +sinogram = radon(image, theta=theta, circle=True) +plt.subplot(122) plt.title("Radon transform\n(Sinogram)") -plt.xlabel("Projection angle (degrees)") -plt.ylabel("Projection axis") -plt.imshow(projections, aspect='auto') - -reconstruction = iradon(projections) -plt.subplot(224) -plt.title("Reconstruction\nfrom sinogram") -plt.imshow(reconstruction, cmap=plt.cm.Greys_r) +plt.xlabel("Projection angle (deg)") +plt.ylabel("Projection position (pixels)") +plt.imshow(sinogram, cmap=plt.cm.Greys_r, + extent=(0, 180, 0, sinogram.shape[0]), aspect='auto') plt.subplots_adjust(hspace=0.4, wspace=0.5) plt.show() + +""" +.. image:: PLOT2RST.current_figure + +Reconstruction with the Filtered Back Projection (FBP) +====================================================== + +The mathematical foundation of the filtered back projection is the Fourier +slice theorem [2]_. It uses Fourier transform of the projection and +interpolation in Fourier space to obtain the 2D Fourier transform of the image, +which is then inverted to form the reconstructed image. The filtered back +projection is among the fastest methods of performing the inverse Radon +transform. The only tunable parameter for the FBP is the filter, which is +applied to the Fourier transformed projections. It may be used to suppress +high frequency noise in the reconstruction. ``skimage`` provides a few +different options for the filter. + +""" + +from skimage.transform import iradon + +reconstruction_fbp = iradon(sinogram, theta=theta, circle=True) +error = reconstruction_fbp - image +print('FBP rms reconstruction error: %.3g' % np.sqrt(np.mean(error**2))) + +imkwargs = dict(vmin=-0.2, vmax=0.2) +plt.figure(figsize=(8, 4.5)) +plt.subplot(121) +plt.title("Reconstruction\nFiltered back projection") +plt.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r) +plt.subplot(122) +plt.title("Reconstruction error\nFiltered back projection") +plt.imshow(reconstruction_fbp - image, cmap=plt.cm.Greys_r, **imkwargs) +plt.show() + +""" +.. image:: PLOT2RST.current_figure + +Reconstruction with the Simultaneous Algebraic Reconstruction Technique +======================================================================= + +Algebraic reconstruction techniques for tomography are based on a +straightforward idea: for a pixelated image the value of a single ray in a +particular projection is simply a sum of all the pixels the ray passes through +on its way through the object. This is a way of expressing the forward Radon +transform. The inverse Radon transform can then be formulated as a (large) set +of linear equations. As each ray passes through a small fraction of the pixels +in the image, this set of equations is sparse, allowing iterative solvers for +sparse linear systems to tackle the system of equations. One iterative method +has been particularly popular, namely Kaczmarz' method [3]_, which has the +property that the solution will approach a least-squares solution of the +equation set. + +The combination of the formulation of the reconstruction problem as a set +of linear equations and an iterative solver makes algebraic techniques +relatively flexible, hence some forms of prior knowledge can be incorporated +with relative ease. + +``skimage`` provides one of the more popular variations of the algebraic +reconstruction techniques: the Simultaneous Algebraic Reconstruction Technique +(SART) [1]_ [4]_. It uses Kaczmarz' method [3]_ as the iterative solver. A good +reconstruction is normally obtained in a single iteration, making the method +computationally effective. Running one or more extra iterations will normally +improve the reconstruction of sharp, high frequency features and reduce the +mean squared error at the expense of increased high frequency noise (the user +will need to decide on what number of iterations is best suited to the problem +at hand. The implementation in ``skimage`` allows prior information of the +form of a lower and upper threshold on the reconstructed values to be supplied +to the reconstruction. + +""" + +from skimage.transform import iradon_sart + +reconstruction_sart = iradon_sart(sinogram, theta=theta) +error = reconstruction_sart - image +print('SART (1 iteration) rms reconstruction error: %.3g' + % np.sqrt(np.mean(error**2))) + +plt.figure(figsize=(8, 8.5)) + +plt.subplot(221) +plt.title("Reconstruction\nSART") +plt.imshow(reconstruction_sart, cmap=plt.cm.Greys_r) +plt.subplot(222) +plt.title("Reconstruction error\nSART") +plt.imshow(reconstruction_sart - image, cmap=plt.cm.Greys_r, **imkwargs) + +# Run a second iteration of SART by supplying the reconstruction +# from the first iteration as an initial estimate +reconstruction_sart2 = iradon_sart(sinogram, theta=theta, + image=reconstruction_sart) +error = reconstruction_sart2 - image +print('SART (2 iterations) rms reconstruction error: %.3g' + % np.sqrt(np.mean(error**2))) + +plt.subplot(223) +plt.title("Reconstruction\nSART, 2 iterations") +plt.imshow(reconstruction_sart2, cmap=plt.cm.Greys_r) +plt.subplot(224) +plt.title("Reconstruction error\nSART, 2 iterations") +plt.imshow(reconstruction_sart2 - image, cmap=plt.cm.Greys_r, **imkwargs) +plt.show() + +""" +.. image:: PLOT2RST.current_figure + + +.. [1] AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging", + IEEE Press 1988. http://www.slaney.org/pct/pct-toc.html +.. [2] Wikipedia, Radon transform, + http://en.wikipedia.org/wiki/Radon_transform#Relationship_with_the_Fourier_transform +.. [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) +.. [4] AH Andersen, AC Kak, "Simultaneous algebraic reconstruction technique + (SART): a superior implementation of the ART algorithm", Ultrasonic + Imaging 6 pp 81--94 (1984) + +""" diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index f079ada4..2cbc7007 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -2,7 +2,7 @@ from ._hough_transform import (hough_circle, hough_ellipse, hough_line, probabilistic_hough_line) from .hough_transform import (hough, probabilistic_hough, hough_peaks, hough_line_peaks) -from .radon_transform import radon, iradon +from .radon_transform import radon, iradon, iradon_sart from .finite_radon_transform import frt2, ifrt2 from .integral import integral_image, integrate from ._geometric import (warp, warp_coords, estimate_transform, @@ -24,6 +24,7 @@ __all__ = ['hough_circle', 'hough_line_peaks', 'radon', 'iradon', + 'iradon_sart', 'frt2', 'ifrt2', 'integral_image', diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx new file mode 100644 index 00000000..91b943b9 --- /dev/null +++ b/skimage/transform/_radon_transform.pyx @@ -0,0 +1,202 @@ +#cython: cdivision=True +#cython: boundscheck=False +#cython: nonecheck=False +#cython: wraparound=False +import numpy as np + +cimport numpy as cnp +cimport cython +from libc.math cimport cos, sin, floor, ceil, sqrt, abs, M_PI + + +cpdef bilinear_ray_sum(cnp.double_t[:, :] image, cnp.double_t theta, + cnp.double_t ray_position): + """ + Compute the projection of an image along a ray. + + Parameters + ---------- + image : 2D array, dtype=float + Image to project. + theta : float + Angle of the projection + ray_position : float + Position of the ray within the projection + + Returns + ------- + projected_value : float + Ray sum along the projection + norm_of_weights : + A measure of how long the ray's path through the reconstruction + circle was + """ + theta = theta / 180. * M_PI + cdef cnp.double_t radius = image.shape[0] // 2 - 1 + cdef cnp.double_t projection_center = image.shape[0] // 2 + cdef cnp.double_t rotation_center = image.shape[0] // 2 + # (s, t) is the (x, y) system rotated by theta + cdef cnp.double_t t = ray_position - projection_center + # s0 is the half-length of the ray's path in the reconstruction circle + cdef cnp.double_t s0 + s0 = sqrt(radius**2 - t**2) if radius**2 >= t**2 else 0. + cdef Py_ssize_t Ns = 2 * ( ceil(2 * s0)) # number of steps + # along the ray + cdef cnp.double_t ray_sum = 0. + cdef cnp.double_t weight_norm = 0. + cdef cnp.double_t ds, dx, dy, x0, y0, x, y, di, dj, + cdef cnp.double_t index_i, index_j, weight + cdef Py_ssize_t k, i, j + if Ns > 0: + # step length between samples + ds = 2 * s0 / Ns + dx = -ds * cos(theta) + dy = -ds * sin(theta) + # point of entry of the ray into the reconstruction circle + x0 = s0 * cos(theta) - t * sin(theta) + y0 = s0 * sin(theta) + t * cos(theta) + for k in range(Ns+1): + x = x0 + k * dx + y = y0 + k * dy + index_i = x + rotation_center + index_j = y + rotation_center + i = floor(index_i) + j = floor(index_j) + di = index_i - floor(index_i) + dj = index_j - floor(index_j) + # Use linear interpolation between values + # Where values fall outside the array, assume zero + if i > 0 and j > 0: + weight = (1. - di) * (1. - dj) * ds + ray_sum += weight * image[i, j] + weight_norm += weight**2 + if i > 0 and j < image.shape[1] - 1: + weight = (1. - di) * dj * ds + ray_sum += weight * image[i, j+1] + weight_norm += weight**2 + if i < image.shape[0] - 1 and j > 0: + weight = di * (1 - dj) * ds + ray_sum += weight * image[i+1, j] + weight_norm += weight**2 + if i < image.shape[0] - 1 and j < image.shape[1] - 1: + weight = di * dj * ds + ray_sum += weight * image[i+1, j+1] + weight_norm += weight**2 + return ray_sum, weight_norm + + +cpdef bilinear_ray_update(cnp.double_t[:, :] image, + cnp.double_t[:, :] image_update, + cnp.double_t theta, cnp.double_t ray_position, + cnp.double_t projected_value): + """ + Compute the update along a ray using bilinear interpolation. + + Parameters + ---------- + image : 2D array, dtype=float + Current reconstruction estimate + image_update : 2D array, dtype=float + Array of same shape as ``image``. Updates will be added to this array. + theta : float + Angle of the projection + ray_position : float + Position of the ray within the projection + projected_value : float + Projected value (from the sinogram) + + Returns + ------- + deviation : + Deviation before updating the image + """ + cdef cnp.double_t ray_sum, weight_norm, deviation + ray_sum, weight_norm = bilinear_ray_sum(image, theta, ray_position) + if weight_norm > 0.: + deviation = -(ray_sum - projected_value) / weight_norm + else: + deviation = 0. + theta = theta / 180. * M_PI + cdef cnp.double_t radius = image.shape[0] // 2 - 1 + cdef cnp.double_t projection_center = image.shape[0] // 2 + cdef cnp.double_t rotation_center = image.shape[0] // 2 + # (s, t) is the (x, y) system rotated by theta + cdef cnp.double_t t = ray_position - projection_center + # s0 is the half-length of the ray's path in the reconstruction circle + cdef cnp.double_t s0 + s0 = sqrt(radius*radius - t*t) if radius**2 >= t**2 else 0. + cdef Py_ssize_t Ns = 2 * ( ceil(2 * s0)) + cdef cnp.double_t hamming_beta = 0.46164 # beta for equiripple Hamming window + + cdef cnp.double_t ds, dx, dy, x0, y0, x, y, di, dj, index_i, index_j + cdef cnp.double_t hamming_window + cdef Py_ssize_t k, i, j + if Ns > 0: + # Step length between samples + ds = 2 * s0 / Ns + dx = -ds * cos(theta) + dy = -ds * sin(theta) + # Point of entry of the ray into the reconstruction circle + x0 = s0 * cos(theta) - t * sin(theta) + y0 = s0 * sin(theta) + t * cos(theta) + for k in range(Ns+1): + x = x0 + k * dx + y = y0 + k * dy + index_i = x + rotation_center + index_j = y + rotation_center + i = floor(index_i) + j = floor(index_j) + di = index_i - floor(index_i) + dj = index_j - floor(index_j) + hamming_window = ((1 - hamming_beta) + - hamming_beta * cos(2 * M_PI * k / (Ns - 1))) + if i > 0 and j > 0: + image_update[i, j] += (deviation * (1. - di) * (1. - dj) + * ds * hamming_window) + if i > 0 and j < image.shape[1] - 1: + image_update[i, j+1] += (deviation * (1. - di) * dj + * ds * hamming_window) + if i < image.shape[0] - 1 and j > 0: + image_update[i+1, j] += (deviation * di * (1 - dj) + * ds * hamming_window) + if i < image.shape[0] - 1 and j < image.shape[1] - 1: + image_update[i+1, j+1] += (deviation * di * dj + * ds * hamming_window) + return deviation + + +@cython.boundscheck(True) +def sart_projection_update(cnp.double_t[:, :] image not None, + cnp.double_t theta, + cnp.double_t[:] projection not None, + cnp.double_t projection_shift=0.): + """ + Compute update to a reconstruction estimate from a single projection + using bilinear interpolation. + + Parameters + ---------- + image : 2D array, dtype=float + Current reconstruction estimate + theta : float + Angle of the projection + projection : 1D array, dtype=float + Projected values, taken from the sinogram + projection_shift : float + Shift the position of the projection by this many pixels before + using it to compute an update to the reconstruction estimate + + Returns + ------- + image_update : 2D array, dtype=float + Array of same shape as ``image`` containing updates that should be + added to ``image`` to improve the reconstruction estimate + """ + cdef cnp.ndarray[cnp.double_t, ndim=2] image_update = np.zeros_like(image) + cdef cnp.double_t ray_position + cdef Py_ssize_t i + for i in range(projection.shape[0]): + ray_position = i + projection_shift + bilinear_ray_update(image, image_update, theta, ray_position, + projection[i]) + return image_update diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index b114adad..944b4a6b 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- """ radon.py - Radon and inverse radon transforms @@ -16,8 +17,9 @@ from __future__ import division import numpy as np from scipy.fftpack import fftshift, fft, ifft from ._warps_cy import _warp_fast +from ._radon_transform import sart_projection_update -__all__ = ["radon", "iradon"] +__all__ = ["radon", "iradon", "iradon_sart"] def radon(image, theta=None, circle=False): @@ -254,3 +256,162 @@ def iradon(radon_image, theta=None, output_size=None, raise ValueError("Unknown interpolation: %s" % interpolation) return reconstructed * np.pi / (2 * len(th)) + + +def order_angles_golden_ratio(theta): + """ + Order angles to reduce the amount of correlated information + in subsequent projections. + + Parameters + ---------- + theta : 1D array of floats + Projection angles in degrees. Duplicate angles are not allowed. + + Returns + ------- + indices : 1D array of unsigned integers + Indices into ``theta`` such that ``theta[indices]`` gives the + approximate golden ratio ordering of the projections. + + Notes + ----- + The method used here is that of the golden ratio introduced + by T. Kohler. + + References: + -Kohler, T. "A projection access scheme for iterative + reconstruction based on the golden section." Nuclear Science + Symposium Conference Record, 2004 IEEE. Vol. 6. IEEE, 2004. + -Winkelmann, Stefanie, et al. "An optimal radial profile order + based on the Golden Ratio for time-resolved MRI." + Medical Imaging, IEEE Transactions on 26.1 (2007): 68-76. + """ + interval = 180 + + def angle_distance(a, b): + difference = a - b + return min(abs(difference % interval), abs(difference % -interval)) + + remaining = list(np.argsort(theta)) # indices into theta + # yield an arbitrary angle to start things off + index = remaining.pop(0) + angle = theta[index] + yield index + # determine subsequent angles using the golden ratio method + angle_increment = interval * (1 - (np.sqrt(5) - 1) / 2) + while remaining: + angle = (angle + angle_increment) % interval + insert_point = np.searchsorted(theta[remaining], angle) + index_below = insert_point - 1 + index_above = 0 if insert_point == len(remaining) else insert_point + distance_below = angle_distance(angle, theta[remaining[index_below]]) + distance_above = angle_distance(angle, theta[remaining[index_above]]) + if distance_below < distance_above: + yield remaining.pop(index_below) + else: + yield remaining.pop(index_above) + + +def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None, + clip=None, relaxation=0.15): + """ + Inverse radon transform + + Reconstruct an image from the radon transform, using a single iteration of + the Simultaneous Algebraic Reconstruction Technique (SART) algorithm. + + Parameters + ---------- + radon_image : 2D array, dtype=float + Image containing radon transform (sinogram). Each column of + the image corresponds to a projection along a different angle. + 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)). + image : 2D array, dtype=float, optional + Image containing an initial reconstruction estimate. Shape of this + array should be ``(radon_image.shape[0], radon_image.shape[0])``. The + default is an array of zeros. + projection_shifts : 1D array, dtype=float + Shift the projections contained in ``radon_image`` (the sinogram) by + this many pixels before reconstructing the image. The i'th value + defines the shift of the i'th column of ``radon_image``. + clip : length-2 sequence of floats + Force all values in the reconstructed tomogram to lie in the range + ``[clip[0], clip[1]]`` + relaxation : float + Relaxation parameter for the update step. A higher value can + improve the convergence rate, but one runs the risk of instabilities. + Values close to or higher than 1 are not recommended. + + Returns + ------- + output : ndarray + Reconstructed image. + + Notes + ----- + Algebraic Reconstruction Techniques are based on formulating the tomography + reconstruction problem as a set of linear equations. Along each ray, + the projected value is the sum of all the values of the cross section along + the ray. A typical feature of SART (and a few other variants of algebraic + techniques) is that it samples the cross section at equidistant points + along the ray, using linear interpolation between the pixel values of the + cross section. The resulting set of linear equations are then solved using + a slightly modified Kaczmarz method. + + When using SART, a single iteration is usually sufficient to obtain a good + reconstruction. Further iterations will tend to enhance high-frequency + information, but will also often increase the noise. + + References: + -AC Kak, M Slaney, "Principles of Computerized Tomographic + Imaging", IEEE Press 1988. + -AH Andersen, AC Kak, "Simultaneous algebraic reconstruction technique + (SART): a superior implementation of the ART algorithm", Ultrasonic + Imaging 6 pp 81--94 (1984) + -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) + -Kohler, T. "A projection access scheme for iterative + reconstruction based on the golden section." Nuclear Science + Symposium Conference Record, 2004 IEEE. Vol. 6. IEEE, 2004. + -Kaczmarz' method, Wikipedia, + http://en.wikipedia.org/wiki/Kaczmarz_method + """ + if radon_image.ndim != 2: + raise ValueError('radon_image must be two dimensional') + reconstructed_shape = (radon_image.shape[0], radon_image.shape[0]) + if theta is None: + theta = np.linspace(0, 180, radon_image.shape[1], endpoint=False) + elif theta.shape != (radon_image.shape[1],): + raise ValueError('Shape of theta (%s) does not match the ' + 'number of projections (%d)' + % (projection_shifts.shape, radon_image.shape[1])) + if image is None: + image = np.zeros(reconstructed_shape, dtype=np.float) + elif image.shape != reconstructed_shape: + raise ValueError('Shape of image (%s) does not match first dimension ' + 'of radon_image (%s)' + % (image.shape, reconstructed_shape)) + if projection_shifts is None: + projection_shifts = np.zeros((radon_image.shape[1],), dtype=np.float) + elif projection_shifts.shape != (radon_image.shape[1],): + raise ValueError('Shape of projection_shifts (%s) does not match the ' + 'number of projections (%d)' + % (projection_shifts.shape, radon_image.shape[1])) + if not clip is None: + if len(clip) != 2: + raise ValueError('clip must be a length-2 sequence') + clip = (float(clip[0]), float(clip[1])) + relaxation = float(relaxation) + + for angle_index in order_angles_golden_ratio(theta): + image_update = sart_projection_update(image, theta[angle_index], + radon_image[:, angle_index], + projection_shifts[angle_index]) + image += relaxation * image_update + if not clip is None: + image = np.clip(image, clip[0], clip[1]) + return image diff --git a/skimage/transform/setup.py b/skimage/transform/setup.py index b0093d87..22f31696 100644 --- a/skimage/transform/setup.py +++ b/skimage/transform/setup.py @@ -15,6 +15,7 @@ def configuration(parent_package='', top_path=None): cython(['_hough_transform.pyx'], working_path=base_path) cython(['_warps_cy.pyx'], working_path=base_path) + cython(['_radon_transform.pyx'], working_path=base_path) config.add_extension('_hough_transform', sources=['_hough_transform.c'], include_dirs=[get_numpy_include_dirs()]) @@ -22,6 +23,10 @@ def configuration(parent_package='', top_path=None): config.add_extension('_warps_cy', sources=['_warps_cy.c'], include_dirs=[get_numpy_include_dirs(), '../_shared']) + config.add_extension('_radon_transform', + sources=['_radon_transform.c'], + include_dirs=[get_numpy_include_dirs()]) + return config if __name__ == '__main__': diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py index f7127aa5..333c298d 100644 --- a/skimage/transform/tests/test_radon_transform.py +++ b/skimage/transform/tests/test_radon_transform.py @@ -1,15 +1,17 @@ -from __future__ import print_function -from __future__ import division +from __future__ import print_function, division import numpy as np from numpy.testing import assert_raises import itertools +import os.path + 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] +__PHANTOM = imread(os.path.join(data_dir, "phantom.png"), + as_grey=True)[::2, ::2] def _get_phantom(): @@ -296,6 +298,85 @@ def test_radon_iradon_circle(): yield check_radon_iradon_circle, interpolation, shape, output_size +def test_order_angles_golden_ratio(): + from skimage.transform.radon_transform import order_angles_golden_ratio + np.random.seed(1231) + lengths = [1, 4, 10, 180] + for l in lengths: + theta_ordered = np.linspace(0, 180, l, endpoint=False) + theta_random = np.random.uniform(0, 180, l) + for theta in (theta_random, theta_ordered): + indices = [x for x in order_angles_golden_ratio(theta)] + # no duplicate indices allowed + assert len(indices) == len(set(indices)) + + +def test_iradon_sart(): + from skimage.io import imread + from skimage import data_dir + from skimage.transform import rescale, radon, iradon_sart + + debug = False + + shepp_logan = imread(os.path.join(data_dir, "phantom.png"), as_grey=True) + image = rescale(shepp_logan, scale=0.4) + theta_ordered = np.linspace(0., 180., image.shape[0], endpoint=False) + theta_missing_wedge = np.linspace(0., 150., image.shape[0], endpoint=True) + for theta, error_factor in ((theta_ordered, 1.), + (theta_missing_wedge, 2.)): + sinogram = radon(image, theta, circle=True) + reconstructed = iradon_sart(sinogram, theta) + + if debug: + from matplotlib import pyplot as plt + plt.figure() + plt.subplot(221) + plt.imshow(image, interpolation='nearest') + plt.subplot(222) + plt.imshow(sinogram, interpolation='nearest') + plt.subplot(223) + plt.imshow(reconstructed, interpolation='nearest') + plt.subplot(224) + plt.imshow(reconstructed - image, interpolation='nearest') + plt.show() + + delta = np.mean(np.abs(reconstructed - image)) + print('delta (1 iteration) =', delta) + assert delta < 0.016 * error_factor + reconstructed = iradon_sart(sinogram, theta, reconstructed) + delta = np.mean(np.abs(reconstructed - image)) + print('delta (2 iterations) =', delta) + assert delta < 0.013 * error_factor + reconstructed = iradon_sart(sinogram, theta, clip=(0, 1)) + delta = np.mean(np.abs(reconstructed - image)) + print('delta (1 iteration, clip) =', delta) + assert delta < 0.015 * error_factor + + np.random.seed(1239867) + shifts = np.random.uniform(-3, 3, sinogram.shape[1]) + x = np.arange(sinogram.shape[0]) + sinogram_shifted = np.vstack(np.interp(x + shifts[i], x, + sinogram[:, i]) + for i in range(sinogram.shape[1])).T + reconstructed = iradon_sart(sinogram_shifted, theta, + projection_shifts=shifts) + if debug: + from matplotlib import pyplot as plt + plt.figure() + plt.subplot(221) + plt.imshow(image, interpolation='nearest') + plt.subplot(222) + plt.imshow(sinogram_shifted, interpolation='nearest') + plt.subplot(223) + plt.imshow(reconstructed, interpolation='nearest') + plt.subplot(224) + plt.imshow(reconstructed - image, interpolation='nearest') + plt.show() + + delta = np.mean(np.abs(reconstructed - image)) + print('delta (1 iteration, shifted sinogram) =', delta) + assert delta < 0.018 * error_factor + if __name__ == "__main__": from numpy.testing import run_module_suite run_module_suite()