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

Merge pull request #546 from jni/3d-slic

Browse Source 
Add support for 3D images in SLIC segmentation
This commit is contained in:
Johannes Schönberger
2013-07-14 23:50:38 -07:00
parent f11d5a1c8d 8736a66504
commit dc27babee2
10 changed files with 445 additions and 132 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/color/__init__.py
+2
View File
@@ -1,4 +1,5 @@
from .colorconv import (convert_colorspace,
guess_spatial_dimensions,
rgb2hsv,
hsv2rgb,
rgb2xyz,
@@ -45,6 +46,7 @@ from .colorlabel import color_dict, label2rgb
__all__ = ['convert_colorspace',
'guess_spatial_dimensions',
'rgb2hsv',
'hsv2rgb',
'rgb2xyz',
skimage/color/colorconv.py
+60 -21
View File
@@ -49,6 +49,37 @@ from ..util import dtype
from skimage._shared.utils import deprecated
def guess_spatial_dimensions(image):
"""Make an educated guess about whether an image has a channels dimension.
Parameters
----------
image : ndarray
The input image.
Returns
-------
spatial_dims : int or None
The number of spatial dimensions of `image`. If ambiguous, the value
is `None`.
Raises
------
ValueError
If the image array has less than two or more than four dimensions.
"""
if image.ndim == 2:
return 2
if image.ndim == 3 and image.shape[-1] != 3:
return 3
if image.ndim == 3 and image.shape[-1] == 3:
return None
if image.ndim == 4 and image.shape[-1] == 3:
return 3
else:
raise ValueError("Expected 2D, 3D, or 4D array, got %iD." % image.ndim)
@deprecated()
def is_rgb(image):
"""Test whether the image is RGB or RGBA.
@@ -72,7 +103,7 @@ def is_gray(image):
Input image.
"""
return np.squeeze(image).ndim == 2
return image.ndim in (2, 3) and not is_rgb(image)
def convert_colorspace(arr, fromspace, tospace):
@@ -129,8 +160,9 @@ def _prepare_colorarray(arr):
"""
arr = np.asanyarray(arr)
if arr.ndim != 3 or arr.shape[2] != 3:
msg = "the input array must be have a shape == (.,.,3))"
if arr.ndim not in [3, 4] or arr.shape[-1] != 3:
msg = ("the input array must be have a shape == (.., ..,[ ..,] 3)), " +
"got (" + (", ".join(map(str, arr.shape))) + ")")
raise ValueError(msg)
return dtype.img_as_float(arr)
@@ -413,12 +445,12 @@ def _convert(matrix, arr):
The converted array.
"""
arr = _prepare_colorarray(arr)
arr = np.swapaxes(arr, 0, 2)
arr = np.swapaxes(arr, 0, -1)
oldshape = arr.shape
arr = np.reshape(arr, (3, -1))
out = np.dot(matrix, arr)
out.shape = oldshape
out = np.swapaxes(out, 2, 0)
out = np.swapaxes(out, -1, 0)
return np.ascontiguousarray(out)
@@ -473,17 +505,19 @@ def rgb2xyz(rgb):
Parameters
----------
rgb : array_like
The image in RGB format, in a 3-D array of shape (.., .., 3).
The image in RGB format, in a 3- or 4-D array of shape
(.., ..,[ ..,] 3).
Returns
-------
out : ndarray
The image in XYZ format, in a 3-D array of shape (.., .., 3).
The image in XYZ format, in a 3- or 4-D array of shape
(.., ..,[ ..,] 3).
Raises
------
ValueError
If `rgb` is not a 3-D array of shape (.., .., 3).
If `rgb` is not a 3- or 4-D array of shape (.., ..,[ ..,] 3).
Notes
-----
@@ -628,23 +662,24 @@ def gray2rgb(image):
Parameters
----------
image : array_like
Input image of shape ``(M, N)``.
Input image of shape ``(M, N [, P])``.
Returns
-------
rgb : ndarray
RGB image of shape ``(M, N, 3)``.
RGB image of shape ``(M, N, [, P], 3)``.
Raises
------
ValueError
If the input is not 2-dimensional.
If the input is not a 2- or 3-dimensional image.
"""
if np.squeeze(image).ndim == 3 and image.shape[2] in (3, 4):
return image
elif image.ndim == 2 or np.squeeze(image).ndim == 2:
return np.dstack((image, image, image))
elif is_gray(image):
image = image[..., np.newaxis]
return np.concatenate((image,)*3, axis=-1)
else:
raise ValueError("Input image expected to be RGB, RGBA or gray.")
@@ -655,17 +690,19 @@ def xyz2lab(xyz):
Parameters
----------
xyz : array_like
The image in XYZ format, in a 3-D array of shape (.., .., 3).
The image in XYZ format, in a 3- or 4-D array of shape
(.., ..,[ ..,] 3).
Returns
-------
out : ndarray
The image in CIE-LAB format, in a 3-D array of shape (.., .., 3).
The image in CIE-LAB format, in a 3- or 4-D array of shape
(.., ..,[ ..,] 3).
Raises
------
ValueError
If `xyz` is not a 3-D array of shape (.., .., 3).
If `xyz` is not a 3-D array of shape (.., ..,[ ..,] 3).
Notes
-----
@@ -695,14 +732,14 @@ def xyz2lab(xyz):
arr[mask] = np.power(arr[mask], 1. / 3.)
arr[~mask] = 7.787 * arr[~mask] + 16. / 116.
x, y, z = arr[:, :, 0], arr[:, :, 1], arr[:, :, 2]
x, y, z = arr[..., 0], arr[..., 1], arr[..., 2]
# Vector scaling
L = (116. * y) - 16.
a = 500.0 * (x - y)
b = 200.0 * (y - z)
return np.dstack([L, a, b])
return np.concatenate([x[..., np.newaxis] for x in [L, a, b]], axis=-1)
def lab2xyz(lab):
@@ -759,17 +796,19 @@ def rgb2lab(rgb):
Parameters
----------
rgb : array_like
The image in RGB format, in a 3-D array of shape (.., .., 3).
The image in RGB format, in a 3- or 4-D array of shape
(.., ..,[ ..,] 3).
Returns
-------
out : ndarray
The image in Lab format, in a 3-D array of shape (.., .., 3).
The image in Lab format, in a 3- or 4-D array of shape
(.., ..,[ ..,] 3).
Raises
------
ValueError
If `rgb` is not a 3-D array of shape (.., .., 3).
If `rgb` is not a 3- or 4-D array of shape (.., ..,[ ..,] 3).
Notes
-----
skimage/color/tests/test_colorconv.py
+15 -1
View File
@@ -33,7 +33,8 @@ from skimage.color import (rgb2hsv, hsv2rgb,
rgb2grey, gray2rgb,
xyz2lab, lab2xyz,
lab2rgb, rgb2lab,
is_rgb, is_gray
is_rgb, is_gray,
guess_spatial_dimensions
)
from skimage import data_dir, data
@@ -41,6 +42,19 @@ from skimage import data_dir, data
import colorsys
def test_guess_spatial_dimensions():
im1 = np.zeros((5, 5))
im2 = np.zeros((5, 5, 5))
im3 = np.zeros((5, 5, 3))
im4 = np.zeros((5, 5, 5, 3))
im5 = np.zeros((5,))
assert_equal(guess_spatial_dimensions(im1), 2)
assert_equal(guess_spatial_dimensions(im2), 3)
assert_equal(guess_spatial_dimensions(im3), None)
assert_equal(guess_spatial_dimensions(im4), 3)
assert_raises(ValueError, guess_spatial_dimensions, im5)
class TestColorconv(TestCase):
img_rgb = imread(os.path.join(data_dir, 'color.png'))
skimage/segmentation/__init__.py
+1 -1
View File
@@ -1,6 +1,6 @@
from .random_walker_segmentation import random_walker
from ._felzenszwalb import felzenszwalb
from ._slic import slic
from .slic_superpixels import slic
from ._quickshift import quickshift
from .boundaries import find_boundaries, visualize_boundaries, mark_boundaries
from ._clear_border import clear_border
skimage/segmentation/_slic.pyx
+64 -104
View File
@@ -2,140 +2,100 @@
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import collections as coll
import numpy as np
from time import time
from scipy import ndimage
cimport numpy as cnp
from ..util import img_as_float
from ..util import img_as_float, regular_grid
from ..color import rgb2lab, gray2rgb
def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
convert2lab=True):
"""Segments image using k-means clustering in Color-(x,y) space.
def _slic_cython(double[:, :, :, ::1] image_zyx,
long[:, :, ::1] nearest_mean,
double[:, :, ::1] distance,
double[:, ::1] means,
float ratio, int max_iter, int n_segments):
"""Helper function for SLIC segmentation.
Parameters
----------
image : (width, height [, 3]) ndarray
Input image.
n_segments : int, optional (default 100)
The (approximate) number of labels in the segmented output image.
ratio: float, optional (default 10)
Balances color-space proximity and image-space proximity.
Higher values give more weight to color-space.
max_iter : int, optional (default 10)
Maximum number of iterations of k-means.
sigma : float, optional (default 1)
Width of Gaussian smoothing kernel for preprocessing. Zero means no
smoothing.
convert2lab : bool, optional (default True)
Whether the input should be converted to Lab colorspace prior to
segmentation. For this purpose, the input is assumed to be RGB. Highly
recommended.
image_zyx : 4D np.ndarray of double, shape (Z, Y, X, 6)
The image with embedded coordinates, that is, `image_zyx[i, j, k]` is
`array([i, j, k, r, g, b])` or `array([i, j, k, L, a, b])`, depending
on the colorspace.
nearest_mean : 3D np.ndarray of long, shape (Z, Y, X)
The (initially empty) label field.
distance : 3D np.ndarray of double, shape (Z, Y, X)
The (initially infinity) array of distances to the nearest centroid.
means : 2D np.ndarray of double, shape (n_segments, 6)
The centroids obtained by SLIC.
ratio : float
The ratio of xyz-space and colorspace in the clustering.
max_iter : int
The maximum number of k-means iterations.
n_segments : int
The approximate/desired number of segments.
Returns
-------
segment_mask : (width, height) ndarray
Integer mask indicating segment labels.
Notes
-----
The image is smoothed using a Gaussian kernel prior to segmentation.
References
----------
.. [1] Radhakrishna Achanta, Appu Shaji, Kevin Smith, Aurelien Lucchi,
Pascal Fua, and Sabine Süsstrunk, SLIC Superpixels Compared to
State-of-the-art Superpixel Methods, TPAMI, May 2012.
Examples
--------
>>> from skimage.segmentation import slic
>>> from skimage.data import lena
>>> img = lena()
>>> segments = slic(img, n_segments=100, ratio=10)
>>> # Increasing the ratio parameter yields more square regions
>>> segments = slic(img, n_segments=100, ratio=20)
nearest_mean : 3D np.ndarray of long, shape (Z, Y, X)
The label field/superpixels found by SLIC.
"""
if image.ndim == 2:
image = gray2rgb(image)
if image.ndim != 3 or image.shape[2] != 3:
ValueError("Only 1- or 3-channel 2D images are supported.")
image = ndimage.gaussian_filter(img_as_float(image), [sigma, sigma, 0])
if convert2lab:
image = rgb2lab(image)
# initialize on grid:
cdef Py_ssize_t height, width
height, width = image.shape[:2]
cdef Py_ssize_t depth, height, width
depth, height, width = (image_zyx.shape[0], image_zyx.shape[1],
image_zyx.shape[2])
# approximate grid size for desired n_segments
cdef Py_ssize_t step = int(np.ceil(np.sqrt(height * width / n_segments)))
grid_y, grid_x = np.mgrid[:height, :width]
means_y = grid_y[::step, ::step]
means_x = grid_x[::step, ::step]
cdef Py_ssize_t step_z, step_y, step_x
slices = regular_grid((depth, height, width), n_segments)
step_z, step_y, step_x = [int(s.step) for s in slices]
means_color = np.zeros((means_y.shape[0], means_y.shape[1], 3))
cdef cnp.ndarray[dtype=cnp.float_t, ndim=2] means \
= np.dstack([means_y, means_x, means_color]).reshape(-1, 5)
cdef cnp.float_t* current_mean
cdef cnp.float_t* mean_entry
n_means = means.shape[0]
# we do the scaling of ratio in the same way as in the SLIC paper
# so the values have the same meaning
ratio = (ratio / float(step)) ** 2
cdef cnp.ndarray[dtype=cnp.float_t, ndim=3] image_yx \
= np.dstack([grid_y, grid_x, image / ratio]).copy("C")
cdef Py_ssize_t i, k, x, y, x_min, x_max, y_min, y_max, changes
cdef Py_ssize_t i, k, x, y, z, x_min, x_max, y_min, y_max, z_min, z_max, \
changes
cdef double dist_mean
cdef cnp.ndarray[dtype=cnp.intp_t, ndim=2] nearest_mean \
= np.zeros((height, width), dtype=np.intp)
cdef cnp.ndarray[dtype=cnp.float_t, ndim=2] distance \
= np.empty((height, width))
cdef cnp.float_t* image_p = <cnp.float_t*> image_yx.data
cdef cnp.float_t* distance_p = <cnp.float_t*> distance.data
cdef cnp.float_t* current_distance
cdef cnp.float_t* current_pixel
cdef double tmp
for i in range(max_iter):
distance.fill(np.inf)
changes = 0
current_mean = <cnp.float_t*> means.data
distance[:, :, :] = np.inf
# assign pixels to means
for k in range(n_means):
# compute windows:
y_min = int(max(current_mean[0] - 2 * step, 0))
y_max = int(min(current_mean[0] + 2 * step, height))
x_min = int(max(current_mean[1] - 2 * step, 0))
x_max = int(min(current_mean[1] + 2 * step, width))
for y in range(y_min, y_max):
current_pixel = &image_p[5 * (y * width + x_min)]
current_distance = &distance_p[y * width + x_min]
for x in range(x_min, x_max):
mean_entry = current_mean
dist_mean = 0
for c in range(5):
# you would think the compiler can optimize the squaring
# itself. mine can't (with O2)
tmp = current_pixel[0] - mean_entry[0]
dist_mean += tmp * tmp
current_pixel += 1
mean_entry += 1
# some precision issue here. Doesnt work if testing ">"
if current_distance[0] - dist_mean > 1e-10:
nearest_mean[y, x] = k
current_distance[0] = dist_mean
changes += 1
current_distance += 1
current_mean += 5
z_min = int(max(means[k, 0] - 2 * step_z, 0))
z_max = int(min(means[k, 0] + 2 * step_z, depth))
y_min = int(max(means[k, 1] - 2 * step_y, 0))
y_max = int(min(means[k, 1] + 2 * step_y, height))
x_min = int(max(means[k, 2] - 2 * step_x, 0))
x_max = int(min(means[k, 2] + 2 * step_x, width))
for z in range(z_min, z_max):
for y in range(y_min, y_max):
for x in range(x_min, x_max):
dist_mean = 0
for c in range(6):
# you would think the compiler can optimize the
# squaring itself. mine can't (with O2)
tmp = image_zyx[z, y, x, c] - means[k, c]
dist_mean += tmp * tmp
# some precision issue here. Doesnt work if testing ">"
if distance[z, y, x] - dist_mean > 1e-10:
nearest_mean[z, y, x] = k
distance[z, y, x] = dist_mean
changes = 1
if changes == 0:
break
# recompute means:
means_list = [np.bincount(nearest_mean.ravel(),
image_yx[:, :, j].ravel()) for j in range(5)]
in_mean = np.bincount(nearest_mean.ravel())
nearest_mean_ravel = np.asarray(nearest_mean).ravel()
means_list = []
for j in range(6):
image_zyx_ravel = np.ascontiguousarray(image_zyx[:, :, :, j]).ravel()
means_list.append(np.bincount(nearest_mean_ravel,
image_zyx_ravel))
in_mean = np.bincount(nearest_mean_ravel)
in_mean[in_mean == 0] = 1
means = (np.vstack(means_list) / in_mean).T.copy("C")
return nearest_mean
return np.ascontiguousarray(nearest_mean)
skimage/segmentation/slic_superpixels.py
+136
View File
@@ -0,0 +1,136 @@
# coding=utf-8
import collections as coll
import numpy as np
from scipy import ndimage
import warnings
from ..util import img_as_float, regular_grid
from ..color import rgb2lab, gray2rgb, guess_spatial_dimensions
from ._slic import _slic_cython
def slic(image, n_segments=100, ratio=10., max_iter=10, sigma=1,
multichannel=None, convert2lab=True):
"""Segments image using k-means clustering in Color-(x,y) space.
Parameters
----------
image : (width, height [, depth] [, 3]) ndarray
Input image, which can be 2D or 3D, and grayscale or multi-channel
(see `multichannel` parameter).
n_segments : int, optional (default: 100)
The (approximate) number of labels in the segmented output image.
ratio: float, optional (default: 10)
Balances color-space proximity and image-space proximity.
Higher values give more weight to color-space.
max_iter : int, optional (default: 10)
Maximum number of iterations of k-means.
sigma : float, optional (default: 1)
Width of Gaussian smoothing kernel for preprocessing. Zero means no
smoothing.
multichannel : bool, optional (default: None)
Whether the last axis of the image is to be interpreted as multiple
channels. Only 3 channels are supported. If `None`, the function will
attempt to guess this, and raise a warning if ambiguous, when the
array has shape (M, N, 3).
convert2lab : bool, optional (default: True)
Whether the input should be converted to Lab colorspace prior to
segmentation. For this purpose, the input is assumed to be RGB. Highly
recommended.
Returns
-------
segment_mask : (width, height) ndarray
Integer mask indicating segment labels.
Raises
------
ValueError
If:
- the image dimension is not 2 or 3 and `multichannel == False`, OR
- the image dimension is not 3 or 4 and `multichannel == True`, OR
- `multichannel == True` and the length of the last dimension of
the image is not 3, OR
Notes
-----
If `sigma > 0` as is default, the image is smoothed using a Gaussian kernel
prior to segmentation.
The image is rescaled to be in [0, 1] prior to processing.
Images of shape (M, N, 3) are interpreted as 2D RGB images by default. To
interpret them as 3D with the last dimension having length 3, use
`multichannel=False`.
References
----------
.. [1] Radhakrishna Achanta, Appu Shaji, Kevin Smith, Aurelien Lucchi,
Pascal Fua, and Sabine Süsstrunk, SLIC Superpixels Compared to
State-of-the-art Superpixel Methods, TPAMI, May 2012.
Examples
--------
>>> from skimage.segmentation import slic
>>> from skimage.data import lena
>>> img = lena()
>>> segments = slic(img, n_segments=100, ratio=10)
>>> # Increasing the ratio parameter yields more square regions
>>> segments = slic(img, n_segments=100, ratio=20)
"""
spatial_dims = guess_spatial_dimensions(image)
if spatial_dims is None and multichannel is None:
msg = ("Images with dimensions (M, N, 3) are interpreted as 2D+RGB" +
" by default. Use `multichannel=False` to interpret as " +
" 3D image with last dimension of length 3.")
warnings.warn(RuntimeWarning(msg))
multichannel = True
elif multichannel is None:
multichannel = (spatial_dims + 1 == image.ndim)
if ((not multichannel and image.ndim not in [2, 3]) or
(multichannel and image.ndim not in [3, 4]) or
(multichannel and image.shape[-1] != 3)):
ValueError("Only 1- or 3-channel 2- or 3-D images are supported.")
image = img_as_float(image)
if not multichannel:
image = gray2rgb(image)
if image.ndim == 3:
# See 2D RGB image as 3D RGB image with Z = 1
image = image[np.newaxis, ...]
if not isinstance(sigma, coll.Iterable):
sigma = np.array([sigma, sigma, sigma, 0])
if (sigma > 0).any():
image = ndimage.gaussian_filter(image, sigma)
if convert2lab:
image = rgb2lab(image)
# initialize on grid:
depth, height, width = image.shape[:3]
# approximate grid size for desired n_segments
grid_z, grid_y, grid_x = np.mgrid[:depth, :height, :width]
slices = regular_grid(image.shape[:3], n_segments)
step_z, step_y, step_x = [int(s.step) for s in slices]
means_z = grid_z[slices]
means_y = grid_y[slices]
means_x = grid_x[slices]
means_color = np.zeros(means_z.shape + (3,))
means = np.concatenate([means_z[..., np.newaxis], means_y[..., np.newaxis],
means_x[..., np.newaxis], means_color
], axis=-1).reshape(-1, 6)
means = np.ascontiguousarray(means)
# we do the scaling of ratio in the same way as in the SLIC paper
# so the values have the same meaning
ratio = (ratio / float(max((step_z, step_y, step_x)))) ** 2
image_zyx = np.concatenate([grid_z[..., np.newaxis],
grid_y[..., np.newaxis],
grid_x[..., np.newaxis],
image / ratio], axis=-1).copy("C")
nearest_mean = np.zeros((depth, height, width), dtype=np.intp)
distance = np.empty((depth, height, width), dtype=np.float)
segment_map = _slic_cython(image_zyx, nearest_mean, distance, means,
ratio, max_iter, n_segments)
if segment_map.shape[0] == 1:
segment_map = segment_map[0]
return segment_map
skimage/segmentation/tests/test_slic.py
+52 -4
View File
@@ -1,9 +1,11 @@
import itertools as it
import warnings
import numpy as np
from numpy.testing import assert_equal, assert_array_equal
from skimage.segmentation import slic
def test_color():
def test_color_2d():
rnd = np.random.RandomState(0)
img = np.zeros((20, 21, 3))
img[:10, :10, 0] = 1
@@ -12,7 +14,9 @@ def test_color():
img += 0.01 * rnd.normal(size=img.shape)
img[img > 1] = 1
img[img < 0] = 0
seg = slic(img, sigma=0, n_segments=4)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
seg = slic(img, sigma=0, n_segments=4)
# we expect 4 segments
assert_equal(len(np.unique(seg)), 4)
@@ -21,7 +25,8 @@ def test_color():
assert_array_equal(seg[:10, 10:], 1)
assert_array_equal(seg[10:, 10:], 3)
def test_gray():
def test_gray_2d():
rnd = np.random.RandomState(0)
img = np.zeros((20, 21))
img[:10, :10] = 0.33
@@ -30,7 +35,7 @@ def test_gray():
img += 0.0033 * rnd.normal(size=img.shape)
img[img > 1] = 1
img[img < 0] = 0
seg = slic(img, sigma=0, n_segments=4, ratio=50.0)
seg = slic(img, sigma=0, n_segments=4, ratio=20.0, multichannel=False)
assert_equal(len(np.unique(seg)), 4)
assert_array_equal(seg[:10, :10], 0)
@@ -38,6 +43,49 @@ def test_gray():
assert_array_equal(seg[:10, 10:], 1)
assert_array_equal(seg[10:, 10:], 3)
def test_color_3d():
rnd = np.random.RandomState(0)
img = np.zeros((20, 21, 22, 3))
slices = []
for dim_size in img.shape[:-1]:
midpoint = dim_size // 2
slices.append((slice(None, midpoint), slice(midpoint, None)))
slices = list(it.product(*slices))
colors = list(it.product(*(([0, 1],) * 3)))
for s, c in zip(slices, colors):
img[s] = c
img += 0.01 * rnd.normal(size=img.shape)
img[img > 1] = 1
img[img < 0] = 0
seg = slic(img, sigma=0, n_segments=8)
assert_equal(len(np.unique(seg)), 8)
for s, c in zip(slices, range(8)):
assert_array_equal(seg[s], c)
def test_gray_3d():
rnd = np.random.RandomState(0)
img = np.zeros((20, 21, 22))
slices = []
for dim_size in img.shape:
midpoint = dim_size // 2
slices.append((slice(None, midpoint), slice(midpoint, None)))
slices = list(it.product(*slices))
shades = np.arange(0, 1.000001, 1.0/7)
for s, sh in zip(slices, shades):
img[s] = sh
img += 0.001 * rnd.normal(size=img.shape)
img[img > 1] = 1
img[img < 0] = 0
seg = slic(img, sigma=0, n_segments=8, ratio=20.0, multichannel=False)
assert_equal(len(np.unique(seg)), 8)
for s, c in zip(slices, range(8)):
assert_array_equal(seg[s], c)
if __name__ == '__main__':
from numpy import testing
testing.run_module_suite()
skimage/util/__init__.py
+3 -1
View File
@@ -11,6 +11,7 @@ if chk < 18: # Use internal version for numpy versions < 1.8.x
else:
from numpy import pad
del numpy, ver, chk
from ._regular_grid import regular_grid
__all__ = ['img_as_float',
@@ -22,4 +23,5 @@ __all__ = ['img_as_float',
'view_as_blocks',
'view_as_windows',
'pad',
'random_noise']
'random_noise',
'regular_grid']
skimage/util/_regular_grid.py
+72
View File
@@ -0,0 +1,72 @@
import numpy as np
def regular_grid(ar_shape, n_points):
"""Find `n_points` regularly spaced along `ar_shape`.
The returned points (as slices) should be as close to cubically-spaced as
possible. Essentially, the points are spaced by the Nth root of the input
array size, where N is the number of dimensions. However, if an array
dimension cannot fit a full step size, it is "discarded", and the
computation is done for only the remaining dimensions.
Parameters
----------
ar_shape : array-like of ints
The shape of the space embedding the grid. `len(ar_shape)` is the
number of dimensions.
n_points : int
The (approximate) number of points to embed in the space.
Returns
-------
slices : list of slice objects
A slice along each dimension of `ar_shape`, such that the intersection
of all the slices give the coordinates of regularly spaced points.
Examples
--------
>>> ar = np.zeros((20, 40))
>>> g = regular_grid(ar.shape, 8)
>>> g
[slice(5.0, None, 10.0), slice(5.0, None, 10.0)]
>>> ar[g] = 1
>>> ar.sum()
8.0
>>> ar = np.zeros((20, 40))
>>> g = regular_grid(ar.shape, 32)
>>> g
[slice(2.0, None, 5.0), slice(2.0, None, 5.0)]
>>> ar[g] = 1
>>> ar.sum()
32.0
>>> ar = np.zeros((3, 20, 40))
>>> g = regular_grid(ar.shape, 8)
>>> g
[slice(1.0, None, 3.0), slice(5.0, None, 10.0), slice(5.0, None, 10.0)]
>>> ar[g] = 1
>>> ar.sum()
8.0
"""
ar_shape = np.asanyarray(ar_shape)
ndim = len(ar_shape)
unsort_dim_idxs = np.argsort(np.argsort(ar_shape))
sorted_dims = np.sort(ar_shape)
space_size = float(np.prod(ar_shape))
if space_size <= n_points:
return [slice(None)] * ndim
stepsizes = (space_size / n_points) ** (1.0 / ndim) * np.ones(ndim)
if (sorted_dims < stepsizes).any():
for dim in range(ndim):
stepsizes[dim] = sorted_dims[dim]
space_size = float(np.prod(sorted_dims[dim+1:]))
stepsizes[dim+1:] = ((space_size / n_points) **
(1.0 / (ndim - dim - 1)))
if (sorted_dims >= stepsizes).all():
break
starts = stepsizes // 2
stepsizes = np.round(stepsizes)
slices = [slice(start, None, step) for
start, step in zip(starts, stepsizes)]
slices = [slices[i] for i in unsort_dim_idxs]
return slices
skimage/util/tests/test_regular_grid.py
+40
View File
@@ -0,0 +1,40 @@
import numpy as np
from numpy.testing import assert_equal
from skimage.util import regular_grid
def test_regular_grid_full():
ar = np.zeros((2, 2))
g = regular_grid(ar, 25)
assert_equal(g, [slice(None, None, None), slice(None, None, None)])
ar[g] = 1
assert_equal(ar.size, ar.sum())
def test_regular_grid_2d_8():
ar = np.zeros((20, 40))
g = regular_grid(ar.shape, 8)
assert_equal(g, [slice(5.0, None, 10.0), slice(5.0, None, 10.0)])
ar[g] = 1
assert_equal(ar.sum(), 8)
def test_regular_grid_2d_32():
ar = np.zeros((20, 40))
g = regular_grid(ar.shape, 32)
assert_equal(g, [slice(2.0, None, 5.0), slice(2.0, None, 5.0)])
ar[g] = 1
assert_equal(ar.sum(), 32)
def test_regular_grid_3d_8():
ar = np.zeros((3, 20, 40))
g = regular_grid(ar.shape, 8)
assert_equal(g, [slice(1.0, None, 3.0), slice(5.0, None, 10.0),
slice(5.0, None, 10.0)])
ar[g] = 1
assert_equal(ar.sum(), 8)
if __name__ == '__main__':
np.testing.run_module_suite()