From 9484afeed17531dbfb747776d05df6c662e5d366 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?= Date: Sun, 2 Jun 2013 21:05:53 +0200 Subject: [PATCH] Add SART tomography reconstruction to radon_transform. --- skimage/transform/_radon_transform.pyx | 170 +++++++++++++++++++++++++ skimage/transform/radon_transform.py | 91 ++++++++++++- skimage/transform/setup.py | 5 + 3 files changed, 265 insertions(+), 1 deletion(-) create mode 100644 skimage/transform/_radon_transform.pyx diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx new file mode 100644 index 00000000..b1bcc300 --- /dev/null +++ b/skimage/transform/_radon_transform.pyx @@ -0,0 +1,170 @@ +#cython: cdivision=True +#cython: boundscheck=True +#cython: nonecheck=True +#cython: wraparound=False +import numpy as np +from numpy import pi + +cimport numpy as cnp +cimport cython +from libc.math cimport cos, sin, floor, ceil, sqrt, abs + + +cpdef bilinear_ray_sum(cnp.ndarray[cnp.double_t, ndim=2] image, double theta, + double ray_position): + '''Compute the projection of an image along a ray. + + Parameters + ---------- + image : 2D array, dtype=float + Image to project. + :param theta: Angle of the projection. + :param ray_position: 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. * pi + cdef double radius = image.shape[0] // 2 - 1 + cdef double projection_center = image.shape[0] // 2 - 1 + cdef double rotation_center = image.shape[0] // 2 + # (s, t) is the (x, y) system rotated by theta + cdef double t = ray_position - projection_center + # s0 is the half-length of the ray's path in the reconstruction circle + cdef double s0 + s0 = sqrt(radius**2 - t**2) if radius**2 >= t**2 else 0. + cdef Py_ssize_t Ns = 2 * int(ceil(2 * s0)) # number of steps along the ray + cdef double ray_sum = 0. + cdef double weight_norm = 0. + cdef double ds, dx, dy, x0, y0, x, y, di, dj, index_i, index_j + 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: + ray_sum += (1. - di) * (1. - dj) * image[i, j] * ds + weight_norm += ((1 - di) * (1 - dj) * ds)**2 + if i > 0 and j < image.shape[1] - 1: + ray_sum += (1. - di) * dj * image[i, j+1] * ds + weight_norm += ((1 - di) * dj * ds)**2 + if i < image.shape[0] - 1 and j > 0: + ray_sum += di * (1 - dj) * image[i+1, j] * ds + weight_norm += (di * (1 - dj) * ds)**2 + if i < image.shape[0] - 1 and j < image.shape[1] - 1: + ray_sum += di * dj * image[i+1, j+1] * ds + weight_norm += (di * dj * ds)**2 + return ray_sum, weight_norm + + +cpdef bilinear_ray_update(cnp.ndarray[cnp.double_t, ndim=2] image, + cnp.ndarray[cnp.double_t, ndim=2] image_update, + double theta, double ray_position, double projected_value): + """Compute the update along a ray using bilinear interpolation. + + Parameters + ---------- + image : + Current reconstruction estimate + image_update : + Array of same shape as ``image``. Updates will be added to this array. + theta : + Angle of the projection + ray_position : + Position of the ray within the projection + projected_value : + Projected value (from the sinogram) + + Returns + ------- + deviation : + Deviation before updating the image + """ + cdef double 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. * pi + cdef double radius = image.shape[0] // 2 - 1 + cdef double projection_center = image.shape[0] // 2 - 1 + cdef double rotation_center = image.shape[0] // 2 + # (s, t) is the (x, y) system rotated by theta + cdef double t = ray_position - projection_center + # s0 is the half-length of the ray's path in the reconstruction circle + cdef double s0 + s0 = sqrt(radius*radius - t*t) if radius**2 >= t**2 else 0. + cdef unsigned int Ns = 2 * int(ceil(2 * s0)) + cdef double hamming_beta = 0.46164 + + cdef double ds, dx, dy, x0, y0, x, y, di, dj, index_i, index_j + cdef double hamming_window + cdef unsigned int 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*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 + + +def sart_projection_update(cnp.ndarray[cnp.double_t, ndim=2] image, \ + double theta, \ + cnp.ndarray[cnp.double_t, ndim=1] projection): + cdef cnp.ndarray[cnp.double_t, ndim=2] image_update = np.zeros_like(image) + cdef unsigned int ray_position + cdef Py_ssize_t i + for i in range(projection.shape[0]): + # TODO: + # ip may differ from i in the future (for alignment of projections) + ray_position = i + 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..55bd1d59 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -16,8 +16,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 +255,91 @@ def iradon(radon_image, theta=None, output_size=None, raise ValueError("Unknown interpolation: %s" % interpolation) return reconstructed * np.pi / (2 * len(th)) + + +def _sart_next_angle(remaining, used): + used = np.array(used) + used.shape = (-1, 1) + remaining = np.array(remaining) + remaining.shape = (1, -1) + time = np.arange(used.shape[0]) + 1 + time.shape = (-1, 1) + tau = 3. + difference = used - remaining + distance = np.minimum(abs(difference % 180), abs(difference % -180)) + #print distance + cost = np.exp(-distance * time / tau).sum(axis=0).squeeze() + next_angle_index = np.argmin(cost) + return remaining[0, next_angle_index] + + +def iradon_sart(radon_image, theta=None, image=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 : 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 (N, M)). + image : array_like, 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. + + 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: + -A. C. Kak, Malcolm 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) + """ + if theta is None: + theta = np.linspace(0, 180, radon_image.shape[1], endpoint=False) + angle_indices = {theta[i]: i for i in range(theta.shape[0])} + reconstructed_shape = (radon_image.shape[0], radon_image.shape[0]) + 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)) + used_angles = [] + while angle_indices: + angle_index = angle_indices.pop(_sart_next_angle(angle_indices.keys(), + used_angles)) + image_update = sart_projection_update(image, theta[angle_index], + radon_image[:, angle_index]) + image += relaxation * image_update + if not clip is None: + image = clip(image, clip[0], clip[1]) + used_angles.append(theta[angle_index]) + + 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__':