wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-09 07:11:31 +08:00
Code Issues Packages Projects Releases Wiki Activity

Do not acquire GIL for iradon_sart

Browse Source
This commit is contained in:
Johannes Schönberger
2015-05-19 18:37:36 -07:00
parent 07f2e4b93f
commit 7a8afbddf3
2 changed files with 84 additions and 76 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/transform/_radon_transform.pyx
+82 -76
View File
@@ -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 * (<Py_ssize_t> ceil(2 * s0)) # number of steps
# along the ray
cdef Py_ssize_t Ns = 2 * (<Py_ssize_t>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 = <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:
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 = <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:
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 * (<Py_ssize_t> ceil(2 * s0))
cdef cnp.double_t hamming_beta = 0.46164 # beta for equiripple Hamming window
cdef Py_ssize_t Ns = 2 * (<Py_ssize_t>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 = <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 * 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 = <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 * 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
skimage/transform/tests/test_radon_transform.py
+2
View File
@@ -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