mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-07 18:05:40 +08:00
Add SART tomography reconstruction to radon_transform.
This commit is contained in:
Jostein Bø Fløystad
@@ -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 = <Py_ssize_t> floor(index_i)
|
||||
j = <Py_ssize_t> 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 = <Py_ssize_t> floor(index_i)
|
||||
j = <Py_ssize_t> 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
|
||||
@@ -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
|
||||
|
||||
@@ -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__':
|
||||
|
||||
Reference in New Issue
Block a user