diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx index 91b943b9..67ddffed 100644 --- a/skimage/transform/_radon_transform.pyx +++ b/skimage/transform/_radon_transform.pyx @@ -40,48 +40,51 @@ cpdef bilinear_ray_sum(cnp.double_t[:, :] image, cnp.double_t theta, # 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 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 + + with nogil: + 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 @@ -89,26 +92,25 @@ 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. + """Compute the update along a ray using bilinear interpolation. Parameters ---------- image : 2D array, dtype=float - Current reconstruction estimate + 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 + Angle of the projection. ray_position : float - Position of the ray within the projection + Position of the ray within the projection. projected_value : float - Projected value (from the sinogram) + Projected value (from the sinogram). Returns ------- deviation : - Deviation before updating the image + 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) @@ -125,43 +127,47 @@ cpdef bilinear_ray_update(cnp.double_t[:, :] image, # 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 Py_ssize_t Ns = 2 * (ceil(2 * s0)) + # beta for equiripple Hamming window + cdef cnp.double_t hamming_beta = 0.46164 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 + + with nogil: + 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 diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py index 9e51ad0a..6df8b7f8 100644 --- a/skimage/transform/tests/test_radon_transform.py +++ b/skimage/transform/tests/test_radon_transform.py @@ -8,6 +8,7 @@ import os.path from skimage.transform import radon, iradon, iradon_sart, rescale from skimage.io import imread from skimage import data_dir +from skimage._shared.testing import test_parallel PHANTOM = imread(os.path.join(data_dir, "phantom.png"), @@ -310,6 +311,7 @@ def test_order_angles_golden_ratio(): assert len(indices) == len(set(indices)) +@test_parallel() def test_iradon_sart(): debug = False