diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx index b1bcc300..e198c4a5 100644 --- a/skimage/transform/_radon_transform.pyx +++ b/skimage/transform/_radon_transform.pyx @@ -157,14 +157,13 @@ cpdef bilinear_ray_update(cnp.ndarray[cnp.double_t, ndim=2] image, def sart_projection_update(cnp.ndarray[cnp.double_t, ndim=2] image, \ double theta, \ - cnp.ndarray[cnp.double_t, ndim=1] projection): + cnp.ndarray[cnp.double_t, ndim=1] projection, + double projection_shift=0.): cdef cnp.ndarray[cnp.double_t, ndim=2] image_update = np.zeros_like(image) - cdef unsigned int ray_position + cdef double 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 + 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 55bd1d59..cc1b76c9 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -273,7 +273,7 @@ def _sart_next_angle(remaining, used): return remaining[0, next_angle_index] -def iradon_sart(radon_image, theta=None, image=None, +def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None, clip=None, relaxation=0.15): """ Inverse radon transform @@ -293,6 +293,10 @@ def iradon_sart(radon_image, theta=None, image=None, 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``. Returns ------- @@ -327,6 +331,8 @@ def iradon_sart(radon_image, theta=None, image=None, reconstructed_shape = (radon_image.shape[0], radon_image.shape[0]) if image is None: image = np.zeros(reconstructed_shape, dtype=np.float) + if projection_shifts is None: + projection_shifts = np.zeros((radon_image.shape[1],), dtype=np.float) elif image.shape != reconstructed_shape: raise ValueError('Shape of image (%s) does not match first dimension ' 'of radon_image (%s)' @@ -336,7 +342,8 @@ def iradon_sart(radon_image, theta=None, image=None, 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]) + radon_image[:, angle_index], + projection_shifts[angle_index]) image += relaxation * image_update if not clip is None: image = clip(image, clip[0], clip[1])