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

Merge pull request #1519 from ahojnnes/nogil

Browse Source 
Release the GIL from Cython
This commit is contained in:
Stefan van der Walt
2015-05-21 10:59:05 -07:00
parent 4e1e861c21 990892fcde
commit bbf3f1c97b
43 changed files with 912 additions and 817 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/_shared/geometry.pxd
+2 -2
View File
@@ -1,6 +1,6 @@
cdef unsigned char point_in_polygon(Py_ssize_t nr_verts, double *xp, double *yp,
double x, double y)
double x, double y) nogil
cdef void points_in_polygon(Py_ssize_t nr_verts, double *xp, double *yp,
Py_ssize_t nr_points, double *x, double *y,
unsigned char *result)
unsigned char *result) nogil
skimage/_shared/geometry.pyx
+3 -2
View File
@@ -5,7 +5,8 @@
cdef inline unsigned char point_in_polygon(Py_ssize_t nr_verts, double *xp,
double *yp, double x, double y):
double *yp, double x,
double y) nogil:
"""Test whether point lies inside a polygon.
Parameters
@@ -33,7 +34,7 @@ cdef inline unsigned char point_in_polygon(Py_ssize_t nr_verts, double *xp,
cdef void points_in_polygon(Py_ssize_t nr_verts, double *xp, double *yp,
Py_ssize_t nr_points, double *x, double *y,
unsigned char *result):
unsigned char *result) nogil:
"""Test whether points lie inside a polygon.
Parameters
skimage/_shared/interpolation.pxd
+11 -10
View File
@@ -5,7 +5,7 @@
from libc.math cimport ceil, floor
cdef inline Py_ssize_t round(double r):
cdef inline Py_ssize_t round(double r) nogil:
return <Py_ssize_t>((r + 0.5) if (r > 0.0) else (r - 0.5))
@@ -13,7 +13,7 @@ cdef inline double nearest_neighbour_interpolation(double* image,
Py_ssize_t rows,
Py_ssize_t cols, double r,
double c, char mode,
double cval):
double cval) nogil:
"""Nearest neighbour interpolation at a given position in the image.
Parameters
@@ -41,7 +41,7 @@ cdef inline double nearest_neighbour_interpolation(double* image,
cdef inline double bilinear_interpolation(double* image, Py_ssize_t rows,
Py_ssize_t cols, double r, double c,
char mode, double cval):
char mode, double cval) nogil:
"""Bilinear interpolation at a given position in the image.
Parameters
@@ -80,7 +80,7 @@ cdef inline double bilinear_interpolation(double* image, Py_ssize_t rows,
return (1 - dr) * top + dr * bottom
cdef inline double quadratic_interpolation(double x, double[3] f):
cdef inline double quadratic_interpolation(double x, double[3] f) nogil:
"""WARNING: Do not use, not implemented correctly.
Quadratic interpolation.
@@ -105,7 +105,8 @@ cdef inline double quadratic_interpolation(double x, double[3] f):
cdef inline double biquadratic_interpolation(double* image, Py_ssize_t rows,
Py_ssize_t cols, double r,
double c, char mode, double cval):
double c, char mode,
double cval) nogil:
"""WARNING: Do not use, not implemented correctly.
Biquadratic interpolation at a given position in the image.
@@ -152,7 +153,7 @@ cdef inline double biquadratic_interpolation(double* image, Py_ssize_t rows,
return quadratic_interpolation(xr, fr)
cdef inline double cubic_interpolation(double x, double[4] f):
cdef inline double cubic_interpolation(double x, double[4] f) nogil:
"""Cubic interpolation.
Parameters
@@ -177,7 +178,7 @@ cdef inline double cubic_interpolation(double x, double[4] f):
cdef inline double bicubic_interpolation(double* image, Py_ssize_t rows,
Py_ssize_t cols, double r, double c,
char mode, double cval):
char mode, double cval) nogil:
"""Bicubic interpolation at a given position in the image.
Interpolation using Catmull-Rom splines, based on the bicubic convolution
@@ -236,7 +237,7 @@ cdef inline double bicubic_interpolation(double* image, Py_ssize_t rows,
cdef inline double get_pixel2d(double* image, Py_ssize_t rows, Py_ssize_t cols,
long r, long c, char mode,
double cval):
double cval) nogil:
"""Get a pixel from the image, taking wrapping mode into consideration.
Parameters
@@ -269,7 +270,7 @@ cdef inline double get_pixel2d(double* image, Py_ssize_t rows, Py_ssize_t cols,
cdef inline double get_pixel3d(double* image, Py_ssize_t rows, Py_ssize_t cols,
Py_ssize_t dims, long r, long c,
long d, char mode, double cval):
long d, char mode, double cval) nogil:
"""Get a pixel from the image, taking wrapping mode into consideration.
Parameters
@@ -302,7 +303,7 @@ cdef inline double get_pixel3d(double* image, Py_ssize_t rows, Py_ssize_t cols,
+ coord_map(dims, d, mode)]
cdef inline Py_ssize_t coord_map(Py_ssize_t dim, long coord, char mode):
cdef inline Py_ssize_t coord_map(Py_ssize_t dim, long coord, char mode) nogil:
"""Wrap a coordinate, according to a given mode.
Parameters
skimage/_shared/transform.pxd
+1 -1
View File
@@ -2,4 +2,4 @@ cimport numpy as cnp
cdef float integrate(float[:, ::1] sat, Py_ssize_t r0, Py_ssize_t c0,
Py_ssize_t r1, Py_ssize_t c1)
Py_ssize_t r1, Py_ssize_t c1) nogil
skimage/_shared/transform.pyx
+1 -1
View File
@@ -6,7 +6,7 @@ cimport numpy as cnp
cdef float integrate(float[:, ::1] sat, Py_ssize_t r0, Py_ssize_t c0,
Py_ssize_t r1, Py_ssize_t c1):
Py_ssize_t r1, Py_ssize_t c1) nogil:
"""
Using a summed area table / integral image, calculate the sum
over a given window.
skimage/draw/_draw.pyx
+34 -32
View File
@@ -6,7 +6,7 @@ import math
import numpy as np
cimport numpy as cnp
from libc.math cimport sqrt, sin, cos, floor, ceil
from libc.math cimport sqrt, sin, cos, floor, ceil, fabs
from .._shared.geometry cimport point_in_polygon
@@ -83,39 +83,41 @@ def line(Py_ssize_t y, Py_ssize_t x, Py_ssize_t y2, Py_ssize_t x2):
cdef Py_ssize_t dy = abs(y2 - y)
cdef Py_ssize_t sx, sy, d, i
if (x2 - x) > 0:
sx = 1
else:
sx = -1
if (y2 - y) > 0:
sy = 1
else:
sy = -1
if dy > dx:
steep = 1
x, y = y, x
dx, dy = dy, dx
sx, sy = sy, sx
d = (2 * dy) - dx
with nogil:
if (x2 - x) > 0:
sx = 1
else:
sx = -1
if (y2 - y) > 0:
sy = 1
else:
sy = -1
if dy > dx:
steep = 1
x, y = y, x
dx, dy = dy, dx
sx, sy = sy, sx
d = (2 * dy) - dx
cdef Py_ssize_t[::1] rr = np.zeros(int(dx) + 1, dtype=np.intp)
cdef Py_ssize_t[::1] cc = np.zeros(int(dx) + 1, dtype=np.intp)
for i in range(dx):
if steep:
rr[i] = x
cc[i] = y
else:
rr[i] = y
cc[i] = x
while d >= 0:
y = y + sy
d = d - (2 * dx)
x = x + sx
d = d + (2 * dy)
with nogil:
for i in range(dx):
if steep:
rr[i] = x
cc[i] = y
else:
rr[i] = y
cc[i] = x
while d >= 0:
y = y + sy
d = d - (2 * dx)
x = x + sx
d = d + (2 * dy)
rr[dx] = y2
cc[dx] = x2
rr[dx] = y2
cc[dx] = x2
return np.asarray(rr), np.asarray(cc)
@@ -187,7 +189,7 @@ def line_aa(Py_ssize_t y1, Py_ssize_t x1, Py_ssize_t y2, Py_ssize_t x2):
while True:
cc.append(x)
rr.append(y)
val.append(1. * abs(err - dx + dy) / <float>(ed))
val.append(1. * fabs(err - dx + dy) / <float>(ed))
e = err
if 2 * e >= -dx:
if x == x2:
@@ -195,7 +197,7 @@ def line_aa(Py_ssize_t y1, Py_ssize_t x1, Py_ssize_t y2, Py_ssize_t x2):
if e + dy < ed:
cc.append(x)
rr.append(y + sign_y)
val.append(1. * abs(e + dy) / <float>(ed))
val.append(1. * fabs(e + dy) / <float>(ed))
err -= dy
x += sign_x
if 2 * e <= dy:
@@ -204,7 +206,7 @@ def line_aa(Py_ssize_t y1, Py_ssize_t x1, Py_ssize_t y2, Py_ssize_t x2):
if dx - e < ed:
cc.append(x)
rr.append(y)
val.append(abs(dx - e) / <float>(ed))
val.append(fabs(dx - e) / <float>(ed))
err += dx
y += sign_y
skimage/draw/tests/test_draw.py
+2
View File
@@ -1,5 +1,6 @@
from numpy.testing import assert_array_equal, assert_equal
import numpy as np
from skimage._shared.testing import test_parallel
from skimage.draw import (set_color, line, line_aa, polygon,
circle, circle_perimeter, circle_perimeter_aa,
@@ -19,6 +20,7 @@ def test_set_color():
assert_array_equal(img, img_)
@test_parallel()
def test_line_horizontal():
img = np.zeros((10, 10))
skimage/feature/_hessian_det_appx.pyx
+23 -21
View File
@@ -6,7 +6,8 @@ import numpy as np
cimport numpy as cnp
cdef inline Py_ssize_t _clip(Py_ssize_t x, Py_ssize_t low, Py_ssize_t high):
cdef inline Py_ssize_t _clip(Py_ssize_t x, Py_ssize_t low,
Py_ssize_t high) nogil:
"""Clips coordinate between high and low.
This method was created so that `hessian_det_appx` does not have to make
@@ -36,7 +37,7 @@ cdef inline Py_ssize_t _clip(Py_ssize_t x, Py_ssize_t low, Py_ssize_t high):
cdef inline cnp.double_t _integ(
cnp.double_t[:, ::1] img, Py_ssize_t r, Py_ssize_t c,
Py_ssize_t rl, Py_ssize_t cl):
Py_ssize_t rl, Py_ssize_t cl) nogil:
"""Integrate over the integral image in the given window
This method was created so that `hessian_det_appx` does not have to make
@@ -123,31 +124,32 @@ def _hessian_matrix_det(cnp.double_t[:, ::1] img, double sigma):
cdef float dxx, dyy, dxy
if not size % 2:
size += 1
with nogil:
if size % 2 == 0:
size += 1
for r in range(height):
for c in range(width):
tl = _integ(img, r - s3, c - s3, s3, s3) # top left
br = _integ(img, r + 1, c + 1, s3, s3) # bottom right
bl = _integ(img, r - s3, c + 1, s3, s3) # bottom left
tr = _integ(img, r + 1, c - s3, s3, s3) # top right
for r in range(height):
for c in range(width):
tl = _integ(img, r - s3, c - s3, s3, s3) # top left
br = _integ(img, r + 1, c + 1, s3, s3) # bottom right
bl = _integ(img, r - s3, c + 1, s3, s3) # bottom left
tr = _integ(img, r + 1, c - s3, s3, s3) # top right
dxy = bl + tr - tl - br
dxy = -dxy * w_i
dxy = bl + tr - tl - br
dxy = -dxy * w_i
mid = _integ(img, r - s3 + 1, c - s2, 2 * s3 - 1, w) # middle box
side = _integ(img, r - s3 + 1, c - s3 / 2, 2 * s3 - 1, s3) # sides
mid = _integ(img, r - s3 + 1, c - s2, 2 * s3 - 1, w) # middle box
side = _integ(img, r - s3 + 1, c - s3 / 2, 2 * s3 - 1, s3) # sides
dxx = mid - 3 * side
dxx = -dxx * w_i
dxx = mid - 3 * side
dxx = -dxx * w_i
mid = _integ(img, r - s2, c - s3 + 1, w, 2 * s3 - 1)
side = _integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1)
mid = _integ(img, r - s2, c - s3 + 1, w, 2 * s3 - 1)
side = _integ(img, r - s3 / 2, c - s3 + 1, s3, 2 * s3 - 1)
dyy = mid - 3 * side
dyy = -dyy * w_i
dyy = mid - 3 * side
dyy = -dyy * w_i
out[r, c] = (dxx * dyy - 0.81 * (dxy * dxy))
out[r, c] = (dxx * dyy - 0.81 * (dxy * dxy))
return out
skimage/feature/_texture.pyx
+146 -135
View File
@@ -8,6 +8,10 @@ from libc.math cimport sin, cos, abs
from .._shared.interpolation cimport bilinear_interpolation, round
cdef extern from "numpy/npy_math.h":
double NAN "NPY_NAN"
def _glcm_loop(cnp.uint8_t[:, ::1] image, double[:] distances,
double[:] angles, Py_ssize_t levels,
cnp.uint32_t[:, :, :, ::1] out):
@@ -36,32 +40,33 @@ def _glcm_loop(cnp.uint8_t[:, ::1] image, double[:] distances,
cnp.uint8_t i, j
cnp.float64_t angle, distance
rows = image.shape[0]
cols = image.shape[1]
with nogil:
rows = image.shape[0]
cols = image.shape[1]
for a_idx in range(len(angles)):
angle = angles[a_idx]
for d_idx in range(len(distances)):
distance = distances[d_idx]
for r in range(rows):
for c in range(cols):
i = image[r, c]
for a_idx in range(angles.shape[0]):
angle = angles[a_idx]
for d_idx in range(distances.shape[0]):
distance = distances[d_idx]
for r in range(rows):
for c in range(cols):
i = image[r, c]
# compute the location of the offset pixel
row = r + <int>round(sin(angle) * distance)
col = c + <int>round(cos(angle) * distance)
# compute the location of the offset pixel
row = r + <int>round(sin(angle) * distance)
col = c + <int>round(cos(angle) * distance)
# make sure the offset is within bounds
if row >= 0 and row < rows and \
col >= 0 and col < cols:
j = image[row, col]
# make sure the offset is within bounds
if row >= 0 and row < rows and \
col >= 0 and col < cols:
j = image[row, col]
if i >= 0 and i < levels and \
j >= 0 and j < levels:
out[i, j, d_idx, a_idx] += 1
if i >= 0 and i < levels and \
j >= 0 and j < levels:
out[i, j, d_idx, a_idx] += 1
cdef inline int _bit_rotate_right(int value, int length):
cdef inline int _bit_rotate_right(int value, int length) nogil:
"""Cyclic bit shift to the right.
Parameters
@@ -132,124 +137,130 @@ def _local_binary_pattern(double[:, ::1] image,
# To compute the variance features
cdef double sum_, var_, texture_i
for r in range(image.shape[0]):
for c in range(image.shape[1]):
for i in range(P):
texture[i] = bilinear_interpolation(&image[0, 0], rows, cols,
r + rp[i], c + cp[i],
'C', 0)
# signed / thresholded texture
for i in range(P):
if texture[i] - image[r, c] >= 0:
signed_texture[i] = 1
else:
signed_texture[i] = 0
lbp = 0
# if method == 'var':
if method == 'V':
# Compute the variance without passing from numpy.
# Following the LBP paper, we're taking a biased estimate
# of the variance (ddof=0)
sum_ = 0.0
var_ = 0.0
with nogil:
for r in range(image.shape[0]):
for c in range(image.shape[1]):
for i in range(P):
texture_i = texture[i]
sum_ += texture_i
var_ += texture_i * texture_i
var_ = (var_ - (sum_ * sum_) / P) / P
if var_ != 0:
lbp = var_
else:
lbp = np.nan
# if method == 'uniform':
elif method == 'U' or method == 'N':
# determine number of 0 - 1 changes
changes = 0
for i in range(P - 1):
changes += abs(signed_texture[i] - signed_texture[i + 1])
if method == 'N':
# Uniform local binary patterns are defined as patterns
# with at most 2 value changes (from 0 to 1 or from 1 to
# 0). Uniform patterns can be caraterized by their number
# `n_ones` of 1. The possible values for `n_ones` range
# from 0 to P.
# Here is an example for P = 4:
# n_ones=0: 0000
# n_ones=1: 0001, 1000, 0100, 0010
# n_ones=2: 0011, 1001, 1100, 0110
# n_ones=3: 0111, 1011, 1101, 1110
# n_ones=4: 1111
#
# For a pattern of size P there are 2 constant patterns
# corresponding to n_ones=0 and n_ones=P. For each other
# value of `n_ones` , i.e n_ones=[1..P-1], there are P
# possible patterns which are related to each other through
# circular permutations. The total number of uniform
# patterns is thus (2 + P * (P - 1)).
# Given any pattern (uniform or not) we must be able to
# associate a unique code:
# 1. Constant patterns patterns (with n_ones=0 and
# n_ones=P) and non uniform patterns are given fixed
# code values.
# 2. Other uniform patterns are indexed considering the
# value of n_ones, and an index called 'rot_index'
# reprenting the number of circular right shifts
# required to obtain the pattern starting from a
# reference position (corresponding to all zeros stacked
# on the right). This number of rotations (or circular
# right shifts) 'rot_index' is efficiently computed by
# considering the positions of the first 1 and the first
# 0 found in the pattern.
if changes <= 2:
# We have a uniform pattern
n_ones = 0 # determines the number of ones
first_one = -1 # position was the first one
first_zero = -1 # position of the first zero
for i in range(P):
if signed_texture[i]:
n_ones += 1
if first_one == -1:
first_one = i
else:
if first_zero == -1:
first_zero = i
if n_ones == 0:
lbp = 0
elif n_ones == P:
lbp = P * (P - 1) + 1
else:
if first_one == 0:
rot_index = n_ones - first_zero
else:
rot_index = P - first_one
lbp = 1 + (n_ones - 1) * P + rot_index
else: # changes > 2
lbp = P * (P - 1) + 2
else: # method != 'N'
if changes <= 2:
for i in range(P):
lbp += signed_texture[i]
texture[i] = bilinear_interpolation(&image[0, 0], rows, cols,
r + rp[i], c + cp[i],
'C', 0)
# signed / thresholded texture
for i in range(P):
if texture[i] - image[r, c] >= 0:
signed_texture[i] = 1
else:
lbp = P + 1
else:
# method == 'default'
for i in range(P):
lbp += signed_texture[i] * weights[i]
signed_texture[i] = 0
# method == 'ror'
if method == 'R':
# shift LBP P times to the right and get minimum value
rotation_chain[0] = <int>lbp
for i in range(1, P):
rotation_chain[i] = \
_bit_rotate_right(rotation_chain[i - 1], P)
lbp = rotation_chain[0]
for i in range(1, P):
lbp = min(lbp, rotation_chain[i])
lbp = 0
output[r, c] = lbp
# if method == 'var':
if method == 'V':
# Compute the variance without passing from numpy.
# Following the LBP paper, we're taking a biased estimate
# of the variance (ddof=0)
sum_ = 0.0
var_ = 0.0
for i in range(P):
texture_i = texture[i]
sum_ += texture_i
var_ += texture_i * texture_i
var_ = (var_ - (sum_ * sum_) / P) / P
if var_ != 0:
lbp = var_
else:
lbp = NAN
# if method == 'uniform':
elif method == 'U' or method == 'N':
# determine number of 0 - 1 changes
changes = 0
for i in range(P - 1):
changes += (signed_texture[i]
- signed_texture[i + 1]) != 0
if method == 'N':
# Uniform local binary patterns are defined as patterns
# with at most 2 value changes (from 0 to 1 or from 1 to
# 0). Uniform patterns can be characterized by their
# number `n_ones` of 1. The possible values for
# `n_ones` range from 0 to P.
#
# Here is an example for P = 4:
# n_ones=0: 0000
# n_ones=1: 0001, 1000, 0100, 0010
# n_ones=2: 0011, 1001, 1100, 0110
# n_ones=3: 0111, 1011, 1101, 1110
# n_ones=4: 1111
#
# For a pattern of size P there are 2 constant patterns
# corresponding to n_ones=0 and n_ones=P. For each other
# value of `n_ones` , i.e n_ones=[1..P-1], there are P
# possible patterns which are related to each other
# through circular permutations. The total number of
# uniform patterns is thus (2 + P * (P - 1)).
# Given any pattern (uniform or not) we must be able to
# associate a unique code:
#
# 1. Constant patterns patterns (with n_ones=0 and
# n_ones=P) and non uniform patterns are given fixed
# code values.
#
# 2. Other uniform patterns are indexed considering the
# value of n_ones, and an index called 'rot_index'
# reprenting the number of circular right shifts
# required to obtain the pattern starting from a
# reference position (corresponding to all zeros stacked
# on the right). This number of rotations (or circular
# right shifts) 'rot_index' is efficiently computed by
# considering the positions of the first 1 and the first
# 0 found in the pattern.
if changes <= 2:
# We have a uniform pattern
n_ones = 0 # determines the number of ones
first_one = -1 # position was the first one
first_zero = -1 # position of the first zero
for i in range(P):
if signed_texture[i]:
n_ones += 1
if first_one == -1:
first_one = i
else:
if first_zero == -1:
first_zero = i
if n_ones == 0:
lbp = 0
elif n_ones == P:
lbp = P * (P - 1) + 1
else:
if first_one == 0:
rot_index = n_ones - first_zero
else:
rot_index = P - first_one
lbp = 1 + (n_ones - 1) * P + rot_index
else: # changes > 2
lbp = P * (P - 1) + 2
else: # method != 'N'
if changes <= 2:
for i in range(P):
lbp += signed_texture[i]
else:
lbp = P + 1
else:
# method == 'default'
for i in range(P):
lbp += signed_texture[i] * weights[i]
# method == 'ror'
if method == 'R':
# shift LBP P times to the right and get minimum value
rotation_chain[0] = <int>lbp
for i in range(1, P):
rotation_chain[i] = \
_bit_rotate_right(rotation_chain[i - 1], P)
lbp = rotation_chain[0]
for i in range(1, P):
lbp = min(lbp, rotation_chain[i])
output[r, c] = lbp
return np.asarray(output)
skimage/feature/brief.py
+6 -5
View File
@@ -140,7 +140,8 @@ class BRIEF(DescriptorExtractor):
"""
assert_nD(image, 2)
np.random.seed(self.sample_seed)
random = np.random.RandomState()
random.seed(self.sample_seed)
image = _prepare_grayscale_input_2D(image)
@@ -151,7 +152,7 @@ class BRIEF(DescriptorExtractor):
desc_size = self.descriptor_size
patch_size = self.patch_size
if self.mode == 'normal':
samples = (patch_size / 5.0) * np.random.randn(desc_size * 8)
samples = (patch_size / 5.0) * random.randn(desc_size * 8)
samples = np.array(samples, dtype=np.int32)
samples = samples[(samples < (patch_size // 2))
& (samples > - (patch_size - 2) // 2)]
@@ -159,9 +160,9 @@ class BRIEF(DescriptorExtractor):
pos1 = samples[:desc_size * 2].reshape(desc_size, 2)
pos2 = samples[desc_size * 2:desc_size * 4].reshape(desc_size, 2)
elif self.mode == 'uniform':
samples = np.random.randint(-(patch_size - 2) // 2,
(patch_size // 2) + 1,
(desc_size * 2, 2))
samples = random.randint(-(patch_size - 2) // 2,
(patch_size // 2) + 1,
(desc_size * 2, 2))
samples = np.array(samples, dtype=np.int32)
pos1, pos2 = np.split(samples, 2)
skimage/feature/brief_cy.pyx
+11 -10
View File
@@ -12,13 +12,14 @@ def _brief_loop(double[:, ::1] image, unsigned char[:, ::1] descriptors,
cdef Py_ssize_t k, d, kr, kc, pr0, pr1, pc0, pc1
for p in range(pos0.shape[0]):
pr0 = pos0[p, 0]
pc0 = pos0[p, 1]
pr1 = pos1[p, 0]
pc1 = pos1[p, 1]
for k in range(keypoints.shape[0]):
kr = keypoints[k, 0]
kc = keypoints[k, 1]
if image[kr + pr0, kc + pc0] < image[kr + pr1, kc + pc1]:
descriptors[k, p] = True
with nogil:
for p in range(pos0.shape[0]):
pr0 = pos0[p, 0]
pc0 = pos0[p, 1]
pr1 = pos1[p, 0]
pc1 = pos1[p, 1]
for k in range(keypoints.shape[0]):
kr = keypoints[k, 0]
kc = keypoints[k, 1]
if image[kr + pr0, kc + pc0] < image[kr + pr1, kc + pc1]:
descriptors[k, p] = True
skimage/feature/censure_cy.pyx
+44 -42
View File
@@ -18,55 +18,57 @@ def _censure_dob_loop(Py_ssize_t n,
cdef Py_ssize_t n2 = 2 * n
cdef double total_weight = inner_weight + outer_weight
# top-left pixel
inner = (integral_img[n2 + n, n2 + n]
+ integral_img[n2 - n - 1, n2 - n - 1]
- integral_img[n2 + n, n2 - n - 1]
- integral_img[n2 - n - 1, n2 + n])
with nogil:
outer = integral_img[2 * n2, 2 * n2]
# top-left pixel
inner = (integral_img[n2 + n, n2 + n]
+ integral_img[n2 - n - 1, n2 - n - 1]
- integral_img[n2 + n, n2 - n - 1]
- integral_img[n2 - n - 1, n2 + n])
filtered_image[n2, n2] = (outer_weight * outer
- total_weight * inner)
outer = integral_img[2 * n2, 2 * n2]
# left column
for i in range(n2 + 1, integral_img.shape[0] - n2):
inner = (integral_img[i + n, n2 + n]
+ integral_img[i - n - 1, n2 - n - 1]
- integral_img[i + n, n2 - n - 1]
- integral_img[i - n - 1, n2 + n])
filtered_image[n2, n2] = (outer_weight * outer
- total_weight * inner)
outer = (integral_img[i + n2, 2 * n2]
- integral_img[i - n2 - 1, 2 * n2])
# left column
for i in range(n2 + 1, integral_img.shape[0] - n2):
inner = (integral_img[i + n, n2 + n]
+ integral_img[i - n - 1, n2 - n - 1]
- integral_img[i + n, n2 - n - 1]
- integral_img[i - n - 1, n2 + n])
filtered_image[i, n2] = (outer_weight * outer
- total_weight * inner)
outer = (integral_img[i + n2, 2 * n2]
- integral_img[i - n2 - 1, 2 * n2])
# top row
for j in range(n2 + 1, integral_img.shape[1] - n2):
inner = (integral_img[n2 + n, j + n]
+ integral_img[n2 - n - 1, j - n - 1]
- integral_img[n2 + n, j - n - 1]
- integral_img[n2 - n - 1, j + n])
filtered_image[i, n2] = (outer_weight * outer
- total_weight * inner)
outer = (integral_img[2 * n2, j + n2]
- integral_img[2 * n2, j - n2 - 1])
filtered_image[n2, j] = (outer_weight * outer
- total_weight * inner)
# remaining block
for i in range(n2 + 1, integral_img.shape[0] - n2):
# top row
for j in range(n2 + 1, integral_img.shape[1] - n2):
inner = (integral_img[i + n, j + n]
+ integral_img[i - n - 1, j - n - 1]
- integral_img[i + n, j - n - 1]
- integral_img[i - n - 1, j + n])
inner = (integral_img[n2 + n, j + n]
+ integral_img[n2 - n - 1, j - n - 1]
- integral_img[n2 + n, j - n - 1]
- integral_img[n2 - n - 1, j + n])
outer = (integral_img[i + n2, j + n2]
+ integral_img[i - n2 - 1, j - n2 - 1]
- integral_img[i + n2, j - n2 - 1]
- integral_img[i - n2 - 1, j + n2])
outer = (integral_img[2 * n2, j + n2]
- integral_img[2 * n2, j - n2 - 1])
filtered_image[i, j] = (outer_weight * outer
- total_weight * inner)
filtered_image[n2, j] = (outer_weight * outer
- total_weight * inner)
# remaining block
for i in range(n2 + 1, integral_img.shape[0] - n2):
for j in range(n2 + 1, integral_img.shape[1] - n2):
inner = (integral_img[i + n, j + n]
+ integral_img[i - n - 1, j - n - 1]
- integral_img[i + n, j - n - 1]
- integral_img[i - n - 1, j + n])
outer = (integral_img[i + n2, j + n2]
+ integral_img[i - n2 - 1, j - n2 - 1]
- integral_img[i + n2, j - n2 - 1]
- integral_img[i - n2 - 1, j + n2])
filtered_image[i, j] = (outer_weight * outer
- total_weight * inner)
skimage/feature/corner_cy.pyx
+71 -66
View File
@@ -5,7 +5,7 @@
import numpy as np
cimport numpy as cnp
from libc.float cimport DBL_MAX
from libc.math cimport atan2
from libc.math cimport atan2, fabs
from ..util import img_as_float, pad
from ..color import rgb2grey
@@ -67,27 +67,30 @@ def corner_moravec(image, Py_ssize_t window_size=1):
cdef double msum, min_msum
cdef Py_ssize_t r, c, br, bc, mr, mc, a, b
for r in range(2 * window_size, rows - 2 * window_size):
for c in range(2 * window_size, cols - 2 * window_size):
min_msum = DBL_MAX
for br in range(r - window_size, r + window_size + 1):
for bc in range(c - window_size, c + window_size + 1):
if br != r and bc != c:
msum = 0
for mr in range(- window_size, window_size + 1):
for mc in range(- window_size, window_size + 1):
msum += (cimage[r + mr, c + mc]
- cimage[br + mr, bc + mc]) ** 2
min_msum = min(msum, min_msum)
out[r, c] = min_msum
with nogil:
for r in range(2 * window_size, rows - 2 * window_size):
for c in range(2 * window_size, cols - 2 * window_size):
min_msum = DBL_MAX
for br in range(r - window_size, r + window_size + 1):
for bc in range(c - window_size, c + window_size + 1):
if br != r and bc != c:
msum = 0
for mr in range(- window_size, window_size + 1):
for mc in range(- window_size, window_size + 1):
msum += (cimage[r + mr, c + mc]
- cimage[br + mr, bc + mc]) ** 2
min_msum = min(msum, min_msum)
out[r, c] = min_msum
return np.asarray(out)
cdef inline double _corner_fast_response(double curr_pixel,
double* circle_intensities,
signed char* bins, signed char state, char n):
signed char* bins, signed char state,
char n) nogil:
cdef char consecutive_count = 0
cdef double curr_response
cdef Py_ssize_t l, m
@@ -97,7 +100,7 @@ cdef inline double _corner_fast_response(double curr_pixel,
if consecutive_count == n:
curr_response = 0
for m in range(16):
curr_response += abs(circle_intensities[m] - curr_pixel)
curr_response += fabs(circle_intensities[m] - curr_pixel)
return curr_response
else:
consecutive_count = 0
@@ -124,49 +127,50 @@ def _corner_fast(double[:, ::1] image, signed char n, double threshold):
cdef double curr_response
for i in range(3, rows - 3):
for j in range(3, cols - 3):
with nogil:
for i in range(3, rows - 3):
for j in range(3, cols - 3):
curr_pixel = image[i, j]
lower_threshold = curr_pixel - threshold
upper_threshold = curr_pixel + threshold
curr_pixel = image[i, j]
lower_threshold = curr_pixel - threshold
upper_threshold = curr_pixel + threshold
for k in range(16):
circle_intensities[k] = image[i + rp[k], j + cp[k]]
if circle_intensities[k] > upper_threshold:
# Brighter pixel
bins[k] = 'b'
elif circle_intensities[k] < lower_threshold:
# Darker pixel
bins[k] = 'd'
else:
# Similar pixel
bins[k] = 's'
for k in range(16):
circle_intensities[k] = image[i + rp[k], j + cp[k]]
if circle_intensities[k] > upper_threshold:
# Brighter pixel
bins[k] = 'b'
elif circle_intensities[k] < lower_threshold:
# Darker pixel
bins[k] = 'd'
else:
# Similar pixel
bins[k] = 's'
# High speed test for n >= 12
if n >= 12:
speed_sum_b = 0
speed_sum_d = 0
for k in range(0, 16, 4):
if bins[k] == 'b':
speed_sum_b += 1
elif bins[k] == 'd':
speed_sum_d += 1
if speed_sum_d < 3 and speed_sum_b < 3:
continue
# High speed test for n >= 12
if n >= 12:
speed_sum_b = 0
speed_sum_d = 0
for k in range(0, 16, 4):
if bins[k] == 'b':
speed_sum_b += 1
elif bins[k] == 'd':
speed_sum_d += 1
if speed_sum_d < 3 and speed_sum_b < 3:
continue
# Test for bright pixels
curr_response = \
_corner_fast_response(curr_pixel, circle_intensities,
bins, 'b', n)
# Test for dark pixels
if curr_response == 0:
# Test for bright pixels
curr_response = \
_corner_fast_response(curr_pixel, circle_intensities,
bins, 'd', n)
bins, 'b', n)
corner_response[i, j] = curr_response
# Test for dark pixels
if curr_response == 0:
curr_response = \
_corner_fast_response(curr_pixel, circle_intensities,
bins, 'd', n)
corner_response[i, j] = curr_response
return np.asarray(corner_response)
@@ -254,22 +258,23 @@ def corner_orientations(image, Py_ssize_t[:, :] corners, mask):
cdef double curr_pixel
cdef double m01, m10, m01_tmp
for i in range(corners.shape[0]):
r0 = corners[i, 0]
c0 = corners[i, 1]
with nogil:
for i in range(corners.shape[0]):
r0 = corners[i, 0]
c0 = corners[i, 1]
m01 = 0
m10 = 0
m01 = 0
m10 = 0
for r in range(mrows):
m01_tmp = 0
for c in range(mcols):
if cmask[r, c]:
curr_pixel = cimage[r0 + r, c0 + c]
m10 += curr_pixel * (c - mcols2)
m01_tmp += curr_pixel
m01 += m01_tmp * (r - mrows2)
for r in range(mrows):
m01_tmp = 0
for c in range(mcols):
if cmask[r, c]:
curr_pixel = cimage[r0 + r, c0 + c]
m10 += curr_pixel * (c - mcols2)
m01_tmp += curr_pixel
m01 += m01_tmp * (r - mrows2)
orientations[i] = atan2(m01, m10)
orientations[i] = atan2(m01, m10)
return np.asarray(orientations)
skimage/feature/orb_cy.pyx
+18 -17
View File
@@ -30,27 +30,28 @@ def _orb_loop(double[:, ::1] image, Py_ssize_t[:, ::1] keypoints,
cdef signed char[:, ::1] cpos0 = POS0
cdef signed char[:, ::1] cpos1 = POS1
for i in range(descriptors.shape[0]):
with nogil:
for i in range(descriptors.shape[0]):
angle = orientations[i]
sin_a = sin(angle)
cos_a = cos(angle)
angle = orientations[i]
sin_a = sin(angle)
cos_a = cos(angle)
kr = keypoints[i, 0]
kc = keypoints[i, 1]
kr = keypoints[i, 0]
kc = keypoints[i, 1]
for j in range(descriptors.shape[1]):
pr0 = cpos0[j, 0]
pc0 = cpos0[j, 1]
pr1 = cpos1[j, 0]
pc1 = cpos1[j, 1]
for j in range(descriptors.shape[1]):
pr0 = cpos0[j, 0]
pc0 = cpos0[j, 1]
pr1 = cpos1[j, 0]
pc1 = cpos1[j, 1]
spr0 = <Py_ssize_t>round(sin_a * pr0 + cos_a * pc0)
spc0 = <Py_ssize_t>round(cos_a * pr0 - sin_a * pc0)
spr1 = <Py_ssize_t>round(sin_a * pr1 + cos_a * pc1)
spc1 = <Py_ssize_t>round(cos_a * pr1 - sin_a * pc1)
spr0 = <Py_ssize_t>round(sin_a * pr0 + cos_a * pc0)
spc0 = <Py_ssize_t>round(cos_a * pr0 - sin_a * pc0)
spr1 = <Py_ssize_t>round(sin_a * pr1 + cos_a * pc1)
spc1 = <Py_ssize_t>round(cos_a * pr1 - sin_a * pc1)
if image[kr + spr0, kc + spc0] < image[kr + spr1, kc + spc1]:
descriptors[i, j] = True
if image[kr + spr0, kc + spc0] < image[kr + spr1, kc + spc1]:
descriptors[i, j] = True
return np.asarray(descriptors)
skimage/feature/tests/test_brief.py
+1
View File
@@ -4,6 +4,7 @@ from skimage import data
from skimage import transform as tf
from skimage.color import rgb2gray
from skimage.feature import BRIEF, corner_peaks, corner_harris
from skimage._shared.testing import test_parallel
def test_color_image_unsupported_error():
skimage/feature/tests/test_censure.py
+2
View File
@@ -2,6 +2,7 @@ import numpy as np
from numpy.testing import assert_array_equal, assert_raises
from skimage.data import moon
from skimage.feature import CENSURE
from skimage._shared.testing import test_parallel
img = moon()
@@ -53,6 +54,7 @@ def test_keypoints_censure_moon_image_dob():
assert_array_equal(expected_scales, detector.scales)
@test_parallel()
def test_keypoints_censure_moon_image_octagon():
"""Verify the actual Censure keypoints and their corresponding scale with
the expected values for Octagon filter."""
skimage/feature/tests/test_corner.py
+5
View File
@@ -6,6 +6,7 @@ from skimage import data
from skimage import img_as_float
from skimage.color import rgb2gray
from skimage.morphology import octagon
from skimage._shared.testing import test_parallel
from skimage.feature import (corner_moravec, corner_harris, corner_shi_tomasi,
corner_subpix, peak_local_max, corner_peaks,
@@ -92,6 +93,7 @@ def test_hessian_matrix_eigvals():
[0, 0, 0, 0, 0]]))
@test_parallel()
def test_hessian_matrix_det():
image = np.zeros((5, 5))
image[2, 2] = 1
@@ -99,6 +101,7 @@ def test_hessian_matrix_det():
assert_almost_equal(det, 0, decimal = 3)
@test_parallel()
def test_square_image():
im = np.zeros((50, 50)).astype(float)
im[:25, :25] = 1.
@@ -280,6 +283,7 @@ def test_corner_fast_image_unsupported_error():
assert_raises(ValueError, corner_fast, img)
@test_parallel()
def test_corner_fast_lena():
img = rgb2gray(data.astronaut())
expected = np.array([[101, 198],
@@ -335,6 +339,7 @@ def test_corner_orientations_even_shape_error():
np.asarray([[7, 7]]), np.ones((4, 4)))
@test_parallel()
def test_corner_orientations_lena():
img = rgb2gray(data.lena())
corners = corner_peaks(corner_fast(img, 11, 0.35))
skimage/feature/tests/test_orb.py
+2
View File
@@ -3,11 +3,13 @@ from numpy.testing import assert_equal, assert_almost_equal, run_module_suite
from skimage.feature import ORB
from skimage import data
from skimage.color import rgb2gray
from skimage._shared.testing import test_parallel
img = rgb2gray(data.lena())
@test_parallel()
def test_keypoints_orb_desired_no_of_keypoints():
detector_extractor = ORB(n_keypoints=10, fast_n=12, fast_threshold=0.20)
detector_extractor.detect(img)
skimage/feature/tests/test_texture.py
+3
View File
@@ -1,5 +1,6 @@
import numpy as np
from skimage.feature import greycomatrix, greycoprops, local_binary_pattern
from skimage._shared.testing import test_parallel
class TestGLCM():
@@ -10,6 +11,7 @@ class TestGLCM():
[0, 2, 2, 2],
[2, 2, 3, 3]], dtype=np.uint8)
@test_parallel()
def test_output_angles(self):
result = greycomatrix(self.image, [1], [0, np.pi / 4, np.pi / 2, 3 * np.pi / 4], 4)
assert result.shape == (4, 4, 1, 4)
@@ -162,6 +164,7 @@ class TestLBP():
[ 0, 255, 30, 34, 255, 24],
[146, 241, 255, 0, 189, 126]], dtype='double')
@test_parallel()
def test_default(self):
lbp = local_binary_pattern(self.image, 8, 1, 'default')
ref = np.array([[ 0, 251, 0, 255, 96, 255],
skimage/filters/_ctmf.pyx
+31 -21
View File
@@ -21,8 +21,8 @@ from libc.string cimport memset
cdef extern from "../_shared/vectorized_ops.h":
void add16(cnp.uint16_t *dest, cnp.uint16_t *src)
void sub16(cnp.uint16_t *dest, cnp.uint16_t *src)
void add16(cnp.uint16_t *dest, cnp.uint16_t *src) nogil
void sub16(cnp.uint16_t *dest, cnp.uint16_t *src) nogil
##############################################################################
@@ -174,14 +174,14 @@ cdef Histograms *allocate_histograms(Py_ssize_t rows,
Py_ssize_t percent,
cnp.uint8_t *data,
cnp.uint8_t *mask,
cnp.uint8_t *output):
cnp.uint8_t *output) nogil:
cdef:
Py_ssize_t adjusted_stripe_length = columns + 2*radius + 1
Py_ssize_t memory_size
void *ptr
Histograms *ph
Py_ssize_t roundoff
Py_ssize_t a
Py_ssize_t a, a_2
SCoord *psc
memory_size = (adjusted_stripe_length *
@@ -292,7 +292,7 @@ cdef Histograms *allocate_histograms(Py_ssize_t rows,
# free_histograms - frees the Histograms structure
#
############################################################################
cdef void free_histograms(Histograms *ph):
cdef void free_histograms(Histograms *ph) nogil:
free(ph.memory)
############################################################################
@@ -301,7 +301,7 @@ cdef void free_histograms(Histograms *ph):
#
############################################################################
cdef void set_stride(Histograms *ph, SCoord *psc):
cdef void set_stride(Histograms *ph, SCoord *psc) nogil:
psc.stride = psc.x * ph.col_stride + psc.y * ph.row_stride
############################################################################
@@ -326,17 +326,23 @@ cdef void set_stride(Histograms *ph, SCoord *psc):
# a column that is "radius" to the left.
#
############################################################################
cdef inline Py_ssize_t tl_br_colidx(Histograms *ph, Py_ssize_t colidx):
@cython.cdivision(True)
cdef inline Py_ssize_t tl_br_colidx(Histograms *ph, Py_ssize_t colidx) nogil:
return (colidx + 3*ph.radius + ph.current_row) % ph.stripe_length
cdef inline Py_ssize_t tr_bl_colidx(Histograms *ph, Py_ssize_t colidx):
@cython.cdivision(True)
cdef inline Py_ssize_t tr_bl_colidx(Histograms *ph, Py_ssize_t colidx) nogil:
return (colidx + 3*ph.radius + ph.row_count-ph.current_row) % \
ph.stripe_length
cdef inline Py_ssize_t leading_edge_colidx(Histograms *ph, Py_ssize_t colidx):
@cython.cdivision(True)
cdef inline Py_ssize_t leading_edge_colidx(Histograms *ph,
Py_ssize_t colidx) nogil:
return (colidx + 5*ph.radius) % ph.stripe_length
cdef inline Py_ssize_t trailing_edge_colidx(Histograms *ph, Py_ssize_t colidx):
@cython.cdivision(True)
cdef inline Py_ssize_t trailing_edge_colidx(Histograms *ph,
Py_ssize_t colidx) nogil:
return (colidx + 3*ph.radius - 1) % ph.stripe_length
############################################################################
@@ -348,7 +354,8 @@ cdef inline Py_ssize_t trailing_edge_colidx(Histograms *ph, Py_ssize_t colidx):
# colidx - the index of the column to add
#
############################################################################
cdef inline void accumulate_coarse_histogram(Histograms *ph, Py_ssize_t colidx):
cdef inline void accumulate_coarse_histogram(Histograms *ph,
Py_ssize_t colidx) nogil:
cdef Py_ssize_t offset
offset = tr_bl_colidx(ph, colidx)
@@ -370,7 +377,8 @@ cdef inline void accumulate_coarse_histogram(Histograms *ph, Py_ssize_t colidx):
# for a given column
#
############################################################################
cdef inline void deaccumulate_coarse_histogram(Histograms *ph, Py_ssize_t colidx):
cdef inline void deaccumulate_coarse_histogram(Histograms *ph,
Py_ssize_t colidx) nogil:
cdef Py_ssize_t offset
#
# The trailing diagonals don't appear until here
@@ -401,7 +409,7 @@ cdef inline void deaccumulate_coarse_histogram(Histograms *ph, Py_ssize_t colidx
############################################################################
cdef inline void accumulate_fine_histogram(Histograms *ph,
Py_ssize_t colidx,
Py_ssize_t fineidx):
Py_ssize_t fineidx) nogil:
cdef:
Py_ssize_t fineoffset = fineidx * 16
Py_ssize_t offset
@@ -425,7 +433,7 @@ cdef inline void accumulate_fine_histogram(Histograms *ph,
############################################################################
cdef inline void deaccumulate_fine_histogram(Histograms *ph,
Py_ssize_t colidx,
Py_ssize_t fineidx):
Py_ssize_t fineidx) nogil:
cdef:
Py_ssize_t fineoffset = fineidx * 16
Py_ssize_t offset
@@ -455,7 +463,7 @@ cdef inline void deaccumulate_fine_histogram(Histograms *ph,
#
############################################################################
cdef inline void accumulate(Histograms *ph):
cdef inline void accumulate(Histograms *ph) nogil:
cdef cnp.int32_t accumulator
accumulate_coarse_histogram(ph, ph.current_column)
deaccumulate_coarse_histogram(ph, ph.current_column)
@@ -480,7 +488,7 @@ cdef inline void accumulate(Histograms *ph):
# to choose remains to be done.
############################################################################
cdef inline void update_fine(Histograms *ph, Py_ssize_t fineidx):
cdef inline void update_fine(Histograms *ph, Py_ssize_t fineidx) nogil:
cdef:
Py_ssize_t first_update_column = ph.last_update_column[fineidx]+1
Py_ssize_t update_limit = ph.current_column+1
@@ -507,7 +515,7 @@ cdef inline void update_histogram(Histograms *ph,
HistogramPiece *hist_piece,
pixel_count_t *pixel_count,
SCoord *last_coord,
SCoord *coord):
SCoord *coord) nogil:
cdef:
Py_ssize_t current_column = ph.current_column
Py_ssize_t current_row = ph.current_row
@@ -548,7 +556,7 @@ cdef inline void update_histogram(Histograms *ph,
# update_current_location - update the histograms at the current location
#
############################################################################
cdef inline void update_current_location(Histograms *ph):
cdef inline void update_current_location(Histograms *ph) nogil:
cdef:
Py_ssize_t current_column = ph.current_column
Py_ssize_t radius = ph.radius
@@ -597,7 +605,7 @@ cdef inline void update_current_location(Histograms *ph):
#
############################################################################
cdef inline cnp.uint8_t find_median(Histograms *ph):
cdef inline cnp.uint8_t find_median(Histograms *ph) nogil:
cdef:
Py_ssize_t pixels_below # of pixels below the median
Py_ssize_t i
@@ -652,12 +660,13 @@ cdef int c_median_filter(Py_ssize_t rows,
Py_ssize_t percent,
cnp.uint8_t *data,
cnp.uint8_t *mask,
cnp.uint8_t *output):
cnp.uint8_t *output) nogil:
cdef:
Histograms *ph
Histogram *phistogram
Py_ssize_t row
Py_ssize_t col
Py_ssize_t end_col
Py_ssize_t i
Py_ssize_t top_left_off
Py_ssize_t top_right_off
@@ -706,7 +715,8 @@ cdef int c_median_filter(Py_ssize_t rows,
# Update locations and coarse accumulator for the octagon
# for points before 0
#
for col in range(-radius, 0 if row >= 0 else columns+radius):
end_col = 0 if row >= 0 else columns + radius
for col in range(-radius, end_col):
ph.current_column = col
ph.current_stride = row * row_stride + col * col_stride
update_current_location(ph)
skimage/filters/rank/bilateral_cy.pyx
+3 -3
View File
@@ -14,7 +14,7 @@ cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef Py_ssize_t bilat_pop = 0
@@ -38,7 +38,7 @@ cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef Py_ssize_t bilat_pop = 0
@@ -57,7 +57,7 @@ cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef Py_ssize_t bilat_pop = 0
skimage/filters/rank/core_cy.pxd
+3 -3
View File
@@ -11,13 +11,13 @@ ctypedef fused dtype_t_out:
double_t
cdef dtype_t _max(dtype_t a, dtype_t b)
cdef dtype_t _min(dtype_t a, dtype_t b)
cdef dtype_t _max(dtype_t a, dtype_t b) nogil
cdef dtype_t _min(dtype_t a, dtype_t b) nogil
cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double,
dtype_t, Py_ssize_t, Py_ssize_t, double,
double, Py_ssize_t, Py_ssize_t),
double, Py_ssize_t, Py_ssize_t) nogil,
dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
skimage/filters/rank/core_cy.pyx
+107 -105
View File
@@ -9,29 +9,29 @@ cimport numpy as cnp
from libc.stdlib cimport malloc, free
cdef inline dtype_t _max(dtype_t a, dtype_t b):
cdef inline dtype_t _max(dtype_t a, dtype_t b) nogil:
return a if a >= b else b
cdef inline dtype_t _min(dtype_t a, dtype_t b):
cdef inline dtype_t _min(dtype_t a, dtype_t b) nogil:
return a if a <= b else b
cdef inline void histogram_increment(Py_ssize_t* histo, double* pop,
dtype_t value):
dtype_t value) nogil:
histo[value] += 1
pop[0] += 1
cdef inline void histogram_decrement(Py_ssize_t* histo, double* pop,
dtype_t value):
dtype_t value) nogil:
histo[value] -= 1
pop[0] -= 1
cdef inline char is_in_mask(Py_ssize_t rows, Py_ssize_t cols,
Py_ssize_t r, Py_ssize_t c,
char* mask):
char* mask) nogil:
"""Check whether given coordinate is within image and mask is true."""
if r < 0 or r > rows - 1 or c < 0 or c > cols - 1:
return 0
@@ -44,7 +44,7 @@ cdef inline char is_in_mask(Py_ssize_t rows, Py_ssize_t cols,
cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double,
dtype_t, Py_ssize_t, Py_ssize_t, double,
double, Py_ssize_t, Py_ssize_t),
double, Py_ssize_t, Py_ssize_t) nogil,
dtype_t[:, ::1] image,
char[:, ::1] selem,
char[:, ::1] mask,
@@ -123,123 +123,125 @@ cdef void _core(void kernel(dtype_t_out*, Py_ssize_t, Py_ssize_t*, double,
t = np.vstack((np.zeros((1, selem.shape[1])), selem))
cdef unsigned char[:, :] t_n = (np.diff(t, axis=0) > 0).view(np.uint8)
for r in range(srows):
for c in range(scols):
if t_e[r, c]:
se_e_r[num_se_e] = r - centre_r
se_e_c[num_se_e] = c - centre_c
num_se_e += 1
if t_w[r, c]:
se_w_r[num_se_w] = r - centre_r
se_w_c[num_se_w] = c - centre_c
num_se_w += 1
if t_n[r, c]:
se_n_r[num_se_n] = r - centre_r
se_n_c[num_se_n] = c - centre_c
num_se_n += 1
if t_s[r, c]:
se_s_r[num_se_s] = r - centre_r
se_s_c[num_se_s] = c - centre_c
num_se_s += 1
with nogil:
for r in range(srows):
for c in range(scols):
rr = r - centre_r
cc = c - centre_c
if selem[r, c]:
for r in range(srows):
for c in range(scols):
if t_e[r, c]:
se_e_r[num_se_e] = r - centre_r
se_e_c[num_se_e] = c - centre_c
num_se_e += 1
if t_w[r, c]:
se_w_r[num_se_w] = r - centre_r
se_w_c[num_se_w] = c - centre_c
num_se_w += 1
if t_n[r, c]:
se_n_r[num_se_n] = r - centre_r
se_n_c[num_se_n] = c - centre_c
num_se_n += 1
if t_s[r, c]:
se_s_r[num_se_s] = r - centre_r
se_s_c[num_se_s] = c - centre_c
num_se_s += 1
for r in range(srows):
for c in range(scols):
rr = r - centre_r
cc = c - centre_c
if selem[r, c]:
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_increment(histo, &pop, image[rr, cc])
r = 0
c = 0
kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin,
p0, p1, s0, s1)
# main loop
r = 0
for even_row in range(0, rows, 2):
# ---> west to east
for c in range(1, cols):
for s in range(num_se_e):
rr = r + se_e_r[s]
cc = c + se_e_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_increment(histo, &pop, image[rr, cc])
for s in range(num_se_w):
rr = r + se_w_r[s]
cc = c + se_w_c[s] - 1
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin,
mid_bin, p0, p1, s0, s1)
r += 1 # pass to the next row
if r >= rows:
break
# ---> north to south
for s in range(num_se_s):
rr = r + se_s_r[s]
cc = c + se_s_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_increment(histo, &pop, image[rr, cc])
r = 0
c = 0
kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin, mid_bin,
p0, p1, s0, s1)
# main loop
r = 0
for even_row in range(0, rows, 2):
# ---> west to east
for c in range(1, cols):
for s in range(num_se_e):
rr = r + se_e_r[s]
cc = c + se_e_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_increment(histo, &pop, image[rr, cc])
for s in range(num_se_w):
rr = r + se_w_r[s]
cc = c + se_w_c[s] - 1
for s in range(num_se_n):
rr = r + se_n_r[s] - 1
cc = c + se_n_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin,
mid_bin, p0, p1, s0, s1)
r += 1 # pass to the next row
if r >= rows:
break
# ---> east to west
for c in range(cols - 2, -1, -1):
for s in range(num_se_w):
rr = r + se_w_r[s]
cc = c + se_w_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_increment(histo, &pop, image[rr, cc])
# ---> north to south
for s in range(num_se_s):
rr = r + se_s_r[s]
cc = c + se_s_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_increment(histo, &pop, image[rr, cc])
for s in range(num_se_e):
rr = r + se_e_r[s]
cc = c + se_e_c[s] + 1
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
for s in range(num_se_n):
rr = r + se_n_r[s] - 1
cc = c + se_n_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin,
mid_bin, p0, p1, s0, s1)
kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin,
mid_bin, p0, p1, s0, s1)
r += 1 # pass to the next row
if r >= rows:
break
# ---> east to west
for c in range(cols - 2, -1, -1):
for s in range(num_se_w):
rr = r + se_w_r[s]
cc = c + se_w_c[s]
# ---> north to south
for s in range(num_se_s):
rr = r + se_s_r[s]
cc = c + se_s_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_increment(histo, &pop, image[rr, cc])
for s in range(num_se_e):
rr = r + se_e_r[s]
cc = c + se_e_c[s] + 1
for s in range(num_se_n):
rr = r + se_n_r[s] - 1
cc = c + se_n_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
kernel(&out[r, c, 0], odepth, histo, pop, image[r, c], max_bin,
mid_bin, p0, p1, s0, s1)
kernel(&out[r, c, 0], odepth, histo, pop, image[r, c],
max_bin, mid_bin, p0, p1, s0, s1)
r += 1 # pass to the next row
if r >= rows:
break
# ---> north to south
for s in range(num_se_s):
rr = r + se_s_r[s]
cc = c + se_s_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_increment(histo, &pop, image[rr, cc])
for s in range(num_se_n):
rr = r + se_n_r[s] - 1
cc = c + se_n_c[s]
if is_in_mask(rows, cols, rr, cc, mask_data):
histogram_decrement(histo, &pop, image[rr, cc])
kernel(&out[r, c, 0], odepth, histo, pop, image[r, c],
max_bin, mid_bin, p0, p1, s0, s1)
# release memory allocated by malloc
free(se_e_r)
free(se_e_c)
free(se_w_r)
free(se_w_c)
free(se_n_r)
free(se_n_c)
free(se_s_r)
free(se_s_c)
free(histo)
# release memory allocated by malloc
free(se_e_r)
free(se_e_c)
free(se_w_r)
free(se_w_c)
free(se_n_r)
free(se_n_c)
free(se_s_r)
free(se_s_c)
free(histo)
skimage/filters/rank/generic_cy.pyx
+19 -19
View File
@@ -14,7 +14,7 @@ cdef inline void _kernel_autolevel(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i, imin, imax, delta
@@ -41,7 +41,7 @@ cdef inline void _kernel_bottomhat(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
@@ -59,7 +59,7 @@ cdef inline void _kernel_equalize(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef Py_ssize_t sum = 0
@@ -79,7 +79,7 @@ cdef inline void _kernel_gradient(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i, imin, imax
@@ -102,7 +102,7 @@ cdef inline void _kernel_maximum(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
@@ -120,7 +120,7 @@ cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef Py_ssize_t mean = 0
@@ -138,7 +138,7 @@ cdef inline void _kernel_subtract_mean(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef Py_ssize_t mean = 0
@@ -156,7 +156,7 @@ cdef inline void _kernel_median(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef double sum = pop / 2.0
@@ -177,7 +177,7 @@ cdef inline void _kernel_minimum(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
@@ -195,7 +195,7 @@ cdef inline void _kernel_modal(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t hmax = 0, imax = 0
@@ -217,7 +217,7 @@ cdef inline void _kernel_enhance_contrast(dtype_t_out* out,
Py_ssize_t max_bin,
Py_ssize_t mid_bin, double p0,
double p1, Py_ssize_t s0,
Py_ssize_t s1):
Py_ssize_t s1) nogil:
cdef Py_ssize_t i, imin, imax
@@ -243,7 +243,7 @@ cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
out[0] = <dtype_t_out>pop
@@ -253,7 +253,7 @@ cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef Py_ssize_t sum = 0
@@ -271,7 +271,7 @@ cdef inline void _kernel_threshold(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef Py_ssize_t mean = 0
@@ -289,7 +289,7 @@ cdef inline void _kernel_tophat(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
@@ -307,7 +307,7 @@ cdef inline void _kernel_noise_filter(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef Py_ssize_t min_i
@@ -334,7 +334,7 @@ cdef inline void _kernel_entropy(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef double e, p
@@ -354,7 +354,7 @@ cdef inline void _kernel_otsu(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef Py_ssize_t max_i
cdef double P, mu1, mu2, q1, new_q1, sigma_b, max_sigma_b
@@ -394,7 +394,7 @@ cdef inline void _kernel_win_hist(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef Py_ssize_t max_i
cdef double scale
skimage/filters/rank/percentile_cy.pyx
+9 -9
View File
@@ -12,7 +12,7 @@ cdef inline void _kernel_autolevel(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i, imin, imax, sum, delta
@@ -46,7 +46,7 @@ cdef inline void _kernel_gradient(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i, imin, imax, sum, delta
@@ -75,7 +75,7 @@ cdef inline void _kernel_mean(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i, sum, mean, n
@@ -101,7 +101,7 @@ cdef inline void _kernel_sum(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i, sum, sum_g, n
@@ -128,7 +128,7 @@ cdef inline void _kernel_subtract_mean(dtype_t_out* out, Py_ssize_t odepth,
Py_ssize_t max_bin,
Py_ssize_t mid_bin, double p0,
double p1, Py_ssize_t s0,
Py_ssize_t s1):
Py_ssize_t s1) nogil:
cdef Py_ssize_t i, sum, mean, n
@@ -156,7 +156,7 @@ cdef inline void _kernel_enhance_contrast(dtype_t_out* out,
Py_ssize_t max_bin,
Py_ssize_t mid_bin, double p0,
double p1, Py_ssize_t s0,
Py_ssize_t s1):
Py_ssize_t s1) nogil:
cdef Py_ssize_t i, imin, imax, sum, delta
@@ -191,7 +191,7 @@ cdef inline void _kernel_percentile(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i
cdef Py_ssize_t sum = 0
@@ -216,7 +216,7 @@ cdef inline void _kernel_pop(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef Py_ssize_t i, sum, n
@@ -237,7 +237,7 @@ cdef inline void _kernel_threshold(dtype_t_out* out, Py_ssize_t odepth,
double pop, dtype_t g,
Py_ssize_t max_bin, Py_ssize_t mid_bin,
double p0, double p1,
Py_ssize_t s0, Py_ssize_t s1):
Py_ssize_t s0, Py_ssize_t s1) nogil:
cdef int i
cdef Py_ssize_t sum = 0
skimage/filters/rank/tests/test_rank.py
+4 -1
View File
@@ -8,6 +8,7 @@ from skimage import data, util, morphology
from skimage.morphology import grey, disk
from skimage.filters import rank
from skimage._shared._warnings import expected_warnings
from skimage._shared.testing import test_parallel
def test_all():
@@ -15,7 +16,9 @@ def test_all():
check_all()
@test_parallel()
def check_all():
np.random.seed(0)
image = np.random.rand(25, 25)
selem = morphology.disk(1)
refs = np.load(os.path.join(skimage.data_dir, "rank_filter_tests.npz"))
@@ -489,7 +492,7 @@ def test_entropy():
assert(np.max(rank.entropy(data, selem)) == 12)
# make sure output is of dtype double
with expected_warnings(['Bitdepth of 11']):
with expected_warnings(['Bitdepth of 11']):
out = rank.entropy(data, np.ones((16, 16), dtype=np.uint8))
assert out.dtype == np.double
skimage/graph/heap.pxd
+7 -7
View File
@@ -19,13 +19,13 @@ cdef class BinaryHeap:
cdef REFERENCE_T *_references
cdef REFERENCE_T _popped_ref
cdef void _add_or_remove_level(self, LEVELS_T add_or_remove)
cdef void _update(self)
cdef void _update_one(self, INDEX_T i)
cdef void _remove(self, INDEX_T i)
cdef void _add_or_remove_level(self, LEVELS_T add_or_remove) nogil
cdef void _update(self) nogil
cdef void _update_one(self, INDEX_T i) nogil
cdef void _remove(self, INDEX_T i) nogil
cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference)
cdef VALUE_T pop_fast(self)
cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference) nogil
cdef VALUE_T pop_fast(self) nogil
cdef class FastUpdateBinaryHeap(BinaryHeap):
cdef readonly REFERENCE_T max_reference
@@ -35,4 +35,4 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
cdef VALUE_T value_of_fast(self, REFERENCE_T reference)
cdef INDEX_T push_if_lower_fast(self, VALUE_T value,
REFERENCE_T reference)
REFERENCE_T reference) nogil
skimage/graph/heap.pyx
+12 -10
View File
@@ -43,7 +43,8 @@ cdef extern from "pyport.h":
cdef VALUE_T inf = Py_HUGE_VAL
# this is handy
cdef inline INDEX_T index_min(INDEX_T a, INDEX_T b): return a if a <= b else b
cdef inline INDEX_T index_min(INDEX_T a, INDEX_T b) nogil:
return a if a <= b else b
cdef class BinaryHeap:
@@ -185,7 +186,7 @@ cdef class BinaryHeap:
## C Maintanance methods
cdef void _add_or_remove_level(self, LEVELS_T add_or_remove):
cdef void _add_or_remove_level(self, LEVELS_T add_or_remove) nogil:
# init indexing ints
cdef INDEX_T i, i1, i2, n
@@ -228,7 +229,7 @@ cdef class BinaryHeap:
self._update()
cdef void _update(self):
cdef void _update(self) nogil:
"""Update the full tree from the bottom up.
This should be done after resizing. """
@@ -251,7 +252,7 @@ cdef class BinaryHeap:
values[ii] = values[i+1]
cdef void _update_one(self, INDEX_T i):
cdef void _update_one(self, INDEX_T i) nogil:
"""Update the tree for one value."""
# shorter name for values
@@ -279,7 +280,7 @@ cdef class BinaryHeap:
i = ii-1
cdef void _remove(self, INDEX_T i1):
cdef void _remove(self, INDEX_T i1) nogil:
"""Remove a value from the heap. By index."""
cdef LEVELS_T levels = self.levels
@@ -314,7 +315,7 @@ cdef class BinaryHeap:
## C Public methods
cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference):
cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference) nogil:
"""The c-method for fast pushing.
Returns the index relative to the start of the last level in the heap."""
@@ -338,7 +339,7 @@ cdef class BinaryHeap:
return count
cdef VALUE_T pop_fast(self):
cdef VALUE_T pop_fast(self) nogil:
"""The c-method for fast popping.
Returns the minimum value. The reference is put in self._popped_ref"""
@@ -543,7 +544,7 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
self._crossref[i] = -1
cdef void _remove(self, INDEX_T i1):
cdef void _remove(self, INDEX_T i1) nogil:
""" Remove a value from the heap. By index. """
cdef LEVELS_T levels = self.levels
cdef INDEX_T count = self.count
@@ -581,7 +582,7 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
self._update_one(i2)
cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference):
cdef INDEX_T push_fast(self, VALUE_T value, REFERENCE_T reference) nogil:
"""The c method for fast pushing.
If the reference is already present, will update its value, otherwise
@@ -611,7 +612,8 @@ cdef class FastUpdateBinaryHeap(BinaryHeap):
self._crossref[reference] = ir
return ir
cdef INDEX_T push_if_lower_fast(self, VALUE_T value, REFERENCE_T reference):
cdef INDEX_T push_if_lower_fast(self, VALUE_T value,
REFERENCE_T reference) nogil:
"""If the reference is already present, will update its value ONLY if
the new value is lower than the old one. If the reference is not
present, this append it. If a value was appended, self._pushed is
skimage/graph/tests/test_heap.py
+2
View File
@@ -1,8 +1,10 @@
import time
import random
import skimage.graph.heap as heap
from skimage._shared.testing import test_parallel
@test_parallel()
def test_heap():
_test_heap(100000, True)
_test_heap(100000, False)
skimage/measure/_ccomp.pxd
+3 -3
View File
@@ -4,6 +4,6 @@ cimport numpy as cnp
DTYPE = cnp.intp
ctypedef cnp.intp_t DTYPE_t
cdef DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n)
cdef void set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root)
cdef void join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m)
cdef DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n) nogil
cdef void set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root) nogil
cdef void join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m) nogil
skimage/measure/_ccomp.pyx
+4 -4
View File
@@ -234,7 +234,7 @@ cdef int ravel_index3D(int x, int y, int z, shape_info *shapeinfo):
# discovered and trees begin to surface.
# When we found out that label 5 and 3 are the same, we assign array[5] = 3.
cdef DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n):
cdef DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n) nogil:
"""Find the root of node n.
Given the example above, for any integer from 1 to 9, 1 is always returned
"""
@@ -244,7 +244,7 @@ cdef DTYPE_t find_root(DTYPE_t *forest, DTYPE_t n):
return root
cdef inline void set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root):
cdef inline void set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root) nogil:
"""
Set all nodes on a path to point to new_root.
Given the example above, given n=9, root=6, it would "reconnect" the tree.
@@ -261,7 +261,7 @@ cdef inline void set_root(DTYPE_t *forest, DTYPE_t n, DTYPE_t root):
forest[n] = root
cdef inline void join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m):
cdef inline void join_trees(DTYPE_t *forest, DTYPE_t n, DTYPE_t m) nogil:
"""Join two trees containing nodes n and m.
If we imagine that in the example tree, the root 1 is not known, we
rather have two disjoint trees with roots 2 and 6.
@@ -416,7 +416,7 @@ def label(input, neighbors=None, background=None, return_num=False,
[0 1 0]
[0 0 1]]
>>> from skimage.measure import label
>>> print(label(x, connectivity=1))
>>> print(label(x, connectivity=1))
[[0 1 1]
[2 3 1]
[2 2 4]]
skimage/segmentation/_felzenszwalb_cy.pyx
+36 -28
View File
@@ -42,7 +42,9 @@ def _felzenszwalb_grey(image, double scale=1, sigma=0.8,
if image.ndim != 2:
raise ValueError("This algorithm works only on single-channel 2d"
"images. Got image of shape %s" % str(image.shape))
image = img_as_float(image)
# rescale scale to behave like in reference implementation
scale = float(scale) / 255.
image = scipy.ndimage.gaussian_filter(image, sigma=sigma)
@@ -55,6 +57,7 @@ def _felzenszwalb_grey(image, double scale=1, sigma=0.8,
cdef cnp.ndarray[cnp.float_t, ndim=1] costs = np.hstack([right_cost.ravel(),
down_cost.ravel(), dright_cost.ravel(),
uright_cost.ravel()]).astype(np.float)
# compute edges between pixels:
height, width = image.shape[:2]
cdef cnp.ndarray[cnp.intp_t, ndim=2] segments \
@@ -65,6 +68,7 @@ def _felzenszwalb_grey(image, double scale=1, sigma=0.8,
uright_edges = np.c_[segments[:-1, 1:].ravel(), segments[1:, :-1].ravel()]
cdef cnp.ndarray[cnp.intp_t, ndim=2] edges \
= np.vstack([right_edges, down_edges, dright_edges, uright_edges])
# initialize data structures for segment size
# and inner cost, then start greedy iteration over edges.
edge_queue = np.argsort(costs)
@@ -75,39 +79,43 @@ def _felzenszwalb_grey(image, double scale=1, sigma=0.8,
cdef cnp.float_t *costs_p = <cnp.float_t*>costs.data
cdef cnp.ndarray[cnp.intp_t, ndim=1] segment_size \
= np.ones(width * height, dtype=np.intp)
# inner cost of segments
cdef cnp.ndarray[cnp.float_t, ndim=1] cint = np.zeros(width * height)
cdef cnp.intp_t seg0, seg1, seg_new, e
cdef float cost, inner_cost0, inner_cost1
# set costs_p back one. we increase it before we use it
# since we might continue before that.
costs_p -= 1
for e in range(costs.size):
seg0 = find_root(segments_p, edges_p[0])
seg1 = find_root(segments_p, edges_p[1])
edges_p += 2
costs_p += 1
if seg0 == seg1:
continue
inner_cost0 = cint[seg0] + scale / segment_size[seg0]
inner_cost1 = cint[seg1] + scale / segment_size[seg1]
if costs_p[0] < min(inner_cost0, inner_cost1):
# update size and cost
join_trees(segments_p, seg0, seg1)
seg_new = find_root(segments_p, seg0)
segment_size[seg_new] = segment_size[seg0] + segment_size[seg1]
cint[seg_new] = costs_p[0]
cdef Py_ssize_t num_costs = costs.size
# postprocessing to remove small segments
edges_p = <cnp.intp_t*>edges.data
for e in range(costs.size):
seg0 = find_root(segments_p, edges_p[0])
seg1 = find_root(segments_p, edges_p[1])
edges_p += 2
if seg0 == seg1:
continue
if segment_size[seg0] < min_size or segment_size[seg1] < min_size:
join_trees(segments_p, seg0, seg1)
with nogil:
# set costs_p back one. we increase it before we use it
# since we might continue before that.
costs_p -= 1
for e in range(num_costs):
seg0 = find_root(segments_p, edges_p[0])
seg1 = find_root(segments_p, edges_p[1])
edges_p += 2
costs_p += 1
if seg0 == seg1:
continue
inner_cost0 = cint[seg0] + scale / segment_size[seg0]
inner_cost1 = cint[seg1] + scale / segment_size[seg1]
if costs_p[0] < min(inner_cost0, inner_cost1):
# update size and cost
join_trees(segments_p, seg0, seg1)
seg_new = find_root(segments_p, seg0)
segment_size[seg_new] = segment_size[seg0] + segment_size[seg1]
cint[seg_new] = costs_p[0]
# postprocessing to remove small segments
edges_p = <cnp.intp_t*>edges.data
for e in range(num_costs):
seg0 = find_root(segments_p, edges_p[0])
seg1 = find_root(segments_p, edges_p[1])
edges_p += 2
if seg0 == seg1:
continue
if segment_size[seg0] < min_size or segment_size[seg1] < min_size:
join_trees(segments_p, seg0, seg1)
# unravel the union find tree
flat = segments.ravel()
skimage/segmentation/_quickshift.pyx
+39 -36
View File
@@ -8,6 +8,7 @@ from itertools import product
cimport numpy as cnp
from libc.math cimport exp, sqrt
from libc.float cimport DBL_MAX
from ..util import img_as_float
from ..color import rgb2lab
@@ -99,19 +100,20 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10,
= np.zeros((height, width))
# compute densities
for r in range(height):
for c in range(width):
r_min, r_max = max(r - w, 0), min(r + w + 1, height)
c_min, c_max = max(c - w, 0), min(c + w + 1, width)
for r_ in range(r_min, r_max):
for c_ in range(c_min, c_max):
dist = 0
for channel in range(channels):
dist += (current_pixel_p[channel] -
image_c[r_, c_, channel])**2
dist += (r - r_)**2 + (c - c_)**2
densities[r, c] += exp(-dist / (2 * kernel_size**2))
current_pixel_p += channels
with nogil:
for r in range(height):
for c in range(width):
r_min, r_max = max(r - w, 0), min(r + w + 1, height)
c_min, c_max = max(c - w, 0), min(c + w + 1, width)
for r_ in range(r_min, r_max):
for c_ in range(c_min, c_max):
dist = 0
for channel in range(channels):
dist += (current_pixel_p[channel] -
image_c[r_, c_, channel])**2
dist += (r - r_)**2 + (c - c_)**2
densities[r, c] += exp(-dist / (2 * kernel_size**2))
current_pixel_p += channels
# this will break ties that otherwise would give us headache
densities += random_state.normal(scale=0.00001, size=(height, width))
@@ -123,29 +125,30 @@ def quickshift(image, ratio=1., float kernel_size=5, max_dist=10,
= np.zeros((height, width))
# find nearest node with higher density
current_pixel_p = image_p
for r in range(height):
for c in range(width):
current_density = densities[r, c]
closest = np.inf
r_min, r_max = max(r - w, 0), min(r + w + 1, height)
c_min, c_max = max(c - w, 0), min(c + w + 1, width)
for r_ in range(r_min, r_max):
for c_ in range(c_min, c_max):
if densities[r_, c_] > current_density:
dist = 0
# We compute the distances twice since otherwise
# we get crazy memory overhead
# (width * height * windowsize**2)
for channel in range(channels):
dist += (current_pixel_p[channel] -
image_c[r_, c_, channel])**2
dist += (r - r_)**2 + (c - c_)**2
if dist < closest:
closest = dist
parent[r, c] = r_ * width + c_
dist_parent[r, c] = sqrt(closest)
current_pixel_p += channels
with nogil:
current_pixel_p = image_p
for r in range(height):
for c in range(width):
current_density = densities[r, c]
closest = DBL_MAX
r_min, r_max = max(r - w, 0), min(r + w + 1, height)
c_min, c_max = max(c - w, 0), min(c + w + 1, width)
for r_ in range(r_min, r_max):
for c_ in range(c_min, c_max):
if densities[r_, c_] > current_density:
dist = 0
# We compute the distances twice since otherwise
# we get crazy memory overhead
# (width * height * windowsize**2)
for channel in range(channels):
dist += (current_pixel_p[channel] -
image_c[r_, c_, channel])**2
dist += (r - r_)**2 + (c - c_)**2
if dist < closest:
closest = dist
parent[r, c] = r_ * width + c_
dist_parent[r, c] = sqrt(closest)
current_pixel_p += channels
dist_parent_flat = dist_parent.ravel()
flat = parent.ravel()
skimage/segmentation/_slic.pyx
+116 -115
View File
@@ -101,88 +101,89 @@ def _slic_cython(double[:, :, :, ::1] image_zyx,
# The reference implementation (Achanta et al.) calls this invxywt
cdef double spatial_weight = float(1) / (step ** 2)
for i in range(max_iter):
change = 0
distance[:, :, :] = DBL_MAX
with nogil:
for i in range(max_iter):
change = 0
distance[:, :, :] = DBL_MAX
# assign pixels to segments
for k in range(n_segments):
# assign pixels to segments
for k in range(n_segments):
# segment coordinate centers
cz = segments[k, 0]
cy = segments[k, 1]
cx = segments[k, 2]
# segment coordinate centers
cz = segments[k, 0]
cy = segments[k, 1]
cx = segments[k, 2]
# compute windows
z_min = <Py_ssize_t>max(cz - 2 * step_z, 0)
z_max = <Py_ssize_t>min(cz + 2 * step_z + 1, depth)
y_min = <Py_ssize_t>max(cy - 2 * step_y, 0)
y_max = <Py_ssize_t>min(cy + 2 * step_y + 1, height)
x_min = <Py_ssize_t>max(cx - 2 * step_x, 0)
x_max = <Py_ssize_t>min(cx + 2 * step_x + 1, width)
# compute windows
z_min = <Py_ssize_t>max(cz - 2 * step_z, 0)
z_max = <Py_ssize_t>min(cz + 2 * step_z + 1, depth)
y_min = <Py_ssize_t>max(cy - 2 * step_y, 0)
y_max = <Py_ssize_t>min(cy + 2 * step_y + 1, height)
x_min = <Py_ssize_t>max(cx - 2 * step_x, 0)
x_max = <Py_ssize_t>min(cx + 2 * step_x + 1, width)
for z in range(z_min, z_max):
dz = (sz * (cz - z)) ** 2
for y in range(y_min, y_max):
dy = (sy * (cy - y)) ** 2
for x in range(x_min, x_max):
dist_center = (dz + dy + (sx * (cx - x)) ** 2) * spatial_weight
dist_color = 0
for c in range(3, n_features):
dist_color += (image_zyx[z, y, x, c - 3]
- segments[k, c]) ** 2
if slic_zero:
dist_center += dist_color / max_dist_color[k]
else:
dist_center += dist_color
for z in range(z_min, z_max):
dz = (sz * (cz - z)) ** 2
for y in range(y_min, y_max):
dy = (sy * (cy - y)) ** 2
for x in range(x_min, x_max):
dist_center = (dz + dy + (sx * (cx - x)) ** 2) * spatial_weight
dist_color = 0
for c in range(3, n_features):
dist_color += (image_zyx[z, y, x, c - 3]
- segments[k, c]) ** 2
if slic_zero:
dist_center += dist_color / max_dist_color[k]
else:
dist_center += dist_color
if distance[z, y, x] > dist_center:
nearest_segments[z, y, x] = k
distance[z, y, x] = dist_center
change = 1
if distance[z, y, x] > dist_center:
nearest_segments[z, y, x] = k
distance[z, y, x] = dist_center
change = 1
# stop if no pixel changed its segment
if change == 0:
break
# stop if no pixel changed its segment
if change == 0:
break
# recompute segment centers
# recompute segment centers
# sum features for all segments
n_segment_elems[:] = 0
segments[:, :] = 0
for z in range(depth):
for y in range(height):
for x in range(width):
k = nearest_segments[z, y, x]
n_segment_elems[k] += 1
segments[k, 0] += z
segments[k, 1] += y
segments[k, 2] += x
for c in range(3, n_features):
segments[k, c] += image_zyx[z, y, x, c - 3]
# divide by number of elements per segment to obtain mean
for k in range(n_segments):
for c in range(n_features):
segments[k, c] /= n_segment_elems[k]
# If in SLICO mode, update the color distance maxima
if slic_zero:
# sum features for all segments
n_segment_elems[:] = 0
segments[:, :] = 0
for z in range(depth):
for y in range(height):
for x in range(width):
k = nearest_segments[z, y, x]
dist_color = 0
n_segment_elems[k] += 1
segments[k, 0] += z
segments[k, 1] += y
segments[k, 2] += x
for c in range(3, n_features):
dist_color += (image_zyx[z, y, x, c - 3] -
segments[k, c]) ** 2
segments[k, c] += image_zyx[z, y, x, c - 3]
# The reference implementation seems to only change
# the color if it increases from previous iteration
if max_dist_color[k] < dist_color:
max_dist_color[k] = dist_color
# divide by number of elements per segment to obtain mean
for k in range(n_segments):
for c in range(n_features):
segments[k, c] /= n_segment_elems[k]
# If in SLICO mode, update the color distance maxima
if slic_zero:
for z in range(depth):
for y in range(height):
for x in range(width):
k = nearest_segments[z, y, x]
dist_color = 0
for c in range(3, n_features):
dist_color += (image_zyx[z, y, x, c - 3] -
segments[k, c]) ** 2
# The reference implementation seems to only change
# the color if it increases from previous iteration
if max_dist_color[k] < dist_color:
max_dist_color[k] = dist_color
return np.asarray(nearest_segments)
@@ -237,54 +238,54 @@ def _enforce_label_connectivity_cython(Py_ssize_t[:, :, ::1] segments,
cdef Py_ssize_t[:, ::1] coord_list = np.zeros((max_size, 3), dtype=np.intp)
# loop through all image
for z in range(depth):
for y in range(height):
for x in range(width):
if connected_segments[z, y, x] >= 0:
continue
# find the component size
adjacent = 0
label = segments[z, y, x]
connected_segments[z, y, x] = current_new_label
current_segment_size = 1
bfs_visited = 0
coord_list[bfs_visited, 0] = z
coord_list[bfs_visited, 1] = y
coord_list[bfs_visited, 2] = x
with nogil:
for z in range(depth):
for y in range(height):
for x in range(width):
if connected_segments[z, y, x] >= 0:
continue
# find the component size
adjacent = 0
label = segments[z, y, x]
connected_segments[z, y, x] = current_new_label
current_segment_size = 1
bfs_visited = 0
coord_list[bfs_visited, 0] = z
coord_list[bfs_visited, 1] = y
coord_list[bfs_visited, 2] = x
#perform a breadth first search to find
# the size of the connected component
while bfs_visited < current_segment_size < max_size:
for i in range(6):
zz = coord_list[bfs_visited, 0] + ddz[i]
yy = coord_list[bfs_visited, 1] + ddy[i]
xx = coord_list[bfs_visited, 2] + ddx[i]
if (0 <= xx < width and
0 <= yy < height and
0 <= zz < depth):
if (segments[zz, yy, xx] == label and
connected_segments[zz, yy, xx] == -1):
connected_segments[zz, yy, xx] = \
current_new_label
coord_list[current_segment_size, 0] = zz
coord_list[current_segment_size, 1] = yy
coord_list[current_segment_size, 2] = xx
current_segment_size += 1
if current_segment_size >= max_size:
break
elif (connected_segments[zz, yy, xx] >= 0 and
connected_segments[zz, yy, xx] != current_new_label):
adjacent = connected_segments[zz, yy, xx]
bfs_visited += 1
#perform a breadth first search to find
# the size of the connected component
while bfs_visited < current_segment_size < max_size:
for i in range(6):
zz = coord_list[bfs_visited, 0] + ddz[i]
yy = coord_list[bfs_visited, 1] + ddy[i]
xx = coord_list[bfs_visited, 2] + ddx[i]
if (0 <= xx < width and
0 <= yy < height and
0 <= zz < depth):
if (segments[zz, yy, xx] == label and
connected_segments[zz, yy, xx] == -1):
connected_segments[zz, yy, xx] = \
current_new_label
coord_list[current_segment_size, 0] = zz
coord_list[current_segment_size, 1] = yy
coord_list[current_segment_size, 2] = xx
current_segment_size += 1
if current_segment_size >= max_size:
break
elif (connected_segments[zz, yy, xx] >= 0 and
connected_segments[zz, yy, xx] != current_new_label):
adjacent = connected_segments[zz, yy, xx]
bfs_visited += 1
# change to an adjacent one, like in the original paper
if current_segment_size < min_size:
for i in range(current_segment_size):
connected_segments[coord_list[i, 0],
coord_list[i, 1],
coord_list[i, 2]] = adjacent
else:
current_new_label += 1
# change to an adjacent one, like in the original paper
if current_segment_size < min_size:
for i in range(current_segment_size):
connected_segments[coord_list[i, 0],
coord_list[i, 1],
coord_list[i, 2]] = adjacent
else:
current_new_label += 1
return np.asarray(connected_segments)
skimage/segmentation/tests/test_felzenszwalb.py
+2 -2
View File
@@ -1,11 +1,11 @@
import numpy as np
from numpy.testing import assert_equal, assert_array_equal
from skimage._shared.testing import assert_greater
from skimage._shared.testing import assert_greater, test_parallel
from skimage.segmentation import felzenszwalb
from skimage import data
@test_parallel()
def test_grey():
# very weak tests. This algorithm is pretty unstable.
img = np.zeros((20, 21))
skimage/segmentation/tests/test_quickshift.py
+2 -1
View File
@@ -1,10 +1,11 @@
import numpy as np
from numpy.testing import assert_equal, assert_array_equal
from nose.tools import assert_true
from skimage._shared.testing import assert_greater
from skimage._shared.testing import assert_greater, test_parallel
from skimage.segmentation import quickshift
@test_parallel()
def test_grey():
rnd = np.random.RandomState(0)
img = np.zeros((20, 21))
skimage/segmentation/tests/test_slic.py
+2
View File
@@ -2,8 +2,10 @@ import itertools as it
import numpy as np
from numpy.testing import assert_equal, assert_raises
from skimage.segmentation import slic
from skimage._shared.testing import test_parallel
@test_parallel()
def test_color_2d():
rnd = np.random.RandomState(0)
img = np.zeros((20, 21, 3))
skimage/transform/_hough_transform.pyx
+28 -23
View File
@@ -79,22 +79,24 @@ def _hough_circle(cnp.ndarray img,
num_circle_pixels = circle_x.size
if normalize:
incr = 1.0 / num_circle_pixels
else:
incr = 1
with nogil:
# For each non zero pixel
for p in range(num_pixels):
# Plug the circle at (px, py),
# its coordinates are (tx, ty)
for c in range(num_circle_pixels):
tx = circle_x[c] + x[p]
ty = circle_y[c] + y[p]
if offset:
acc[i, tx, ty] += incr
elif 0 <= tx < xmax and 0 <= ty < ymax:
acc[i, tx, ty] += incr
if normalize:
incr = 1.0 / num_circle_pixels
else:
incr = 1
# For each non zero pixel
for p in range(num_pixels):
# Plug the circle at (px, py),
# its coordinates are (tx, ty)
for c in range(num_circle_pixels):
tx = circle_x[c] + x[p]
ty = circle_y[c] + y[p]
if offset:
acc[i, tx, ty] += incr
elif 0 <= tx < xmax and 0 <= ty < ymax:
acc[i, tx, ty] += incr
return acc
@@ -306,14 +308,17 @@ def hough_line(cnp.ndarray img,
# finally, run the transform
cdef Py_ssize_t nidxs, nthetas, i, j, x, y, accum_idx
nidxs = y_idxs.shape[0] # x and y are the same shape
nthetas = theta.shape[0]
for i in range(nidxs):
x = x_idxs[i]
y = y_idxs[i]
for j in range(nthetas):
accum_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
accum[accum_idx, j] += 1
with nogil:
nidxs = y_idxs.shape[0] # x and y are the same shape
nthetas = theta.shape[0]
for i in range(nidxs):
x = x_idxs[i]
y = y_idxs[i]
for j in range(nthetas):
accum_idx = <int>round((ctheta[j] * x + stheta[j] * y)) + offset
accum[accum_idx, j] += 1
return accum, theta, bins
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/_warps_cy.pyx
+8 -7
View File
@@ -11,7 +11,7 @@ from .._shared.interpolation cimport (nearest_neighbour_interpolation,
cdef inline void _matrix_transform(double x, double y, double* H, double *x_,
double *y_):
double *y_) nogil:
"""Apply a homography to a coordinate.
Parameters
@@ -102,7 +102,7 @@ def _warp_fast(cnp.ndarray image, cnp.ndarray H, output_shape=None,
cdef Py_ssize_t cols = img.shape[1]
cdef double (*interp_func)(double*, Py_ssize_t, Py_ssize_t, double, double,
char, double)
char, double) nogil
if order == 0:
interp_func = nearest_neighbour_interpolation
elif order == 1:
@@ -112,10 +112,11 @@ def _warp_fast(cnp.ndarray image, cnp.ndarray H, output_shape=None,
elif order == 3:
interp_func = bicubic_interpolation
for tfr in range(out_r):
for tfc in range(out_c):
_matrix_transform(tfc, tfr, &M[0, 0], &c, &r)
out[tfr, tfc] = interp_func(&img[0, 0], rows, cols, r, c,
mode_c, cval)
with nogil:
for tfr in range(out_r):
for tfc in range(out_c):
_matrix_transform(tfc, tfr, &M[0, 0], &c, &r)
out[tfr, tfc] = interp_func(&img[0, 0], rows, cols, r, c,
mode_c, cval)
return np.asarray(out)
skimage/transform/tests/test_hough_transform.py
+3
View File
@@ -4,6 +4,7 @@ from numpy.testing import assert_almost_equal, assert_equal
import skimage.transform as tf
from skimage.draw import line, circle_perimeter, ellipse_perimeter
from skimage._shared._warnings import expected_warnings
from skimage._shared.testing import test_parallel
def append_desc(func, description):
@@ -15,6 +16,7 @@ def append_desc(func, description):
return func
@test_parallel()
def test_hough_line():
# Generate a test image
img = np.zeros((100, 150), dtype=int)
@@ -129,6 +131,7 @@ def test_hough_line_peaks_num():
min_angle=0, num_peaks=1)[0]) == 1
@test_parallel()
def test_hough_circle():
# Prepare picture
img = np.zeros((120, 100), dtype=int)
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
skimage/transform/tests/test_warps.py
+2
View File
@@ -11,6 +11,7 @@ from skimage.transform import (warp, warp_coords, rotate, resize, rescale,
from skimage import transform as tf, data, img_as_float
from skimage.color import rgb2gray
from skimage._shared._warnings import expected_warnings
from skimage._shared.testing import test_parallel
np.random.seed(0)
@@ -41,6 +42,7 @@ def test_warp_callable():
assert_almost_equal(outx, refx)
@test_parallel()
def test_warp_matrix():
x = np.zeros((5, 5), dtype=np.double)
x[2, 2] = 1