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

Merge pull request #275 from ahojnnes/warps

Browse Source 
ENH: Refactor image warps.
This commit is contained in:
Stefan van der Walt
2012-08-31 18:12:08 -07:00
parent 06092023b8 b2e4fd6f32
commit c8c1d5c7e7
11 changed files with 487 additions and 138 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/_shared/interpolation.pxd
+17 -6
View File
@@ -1,13 +1,24 @@
cdef inline double nearest_neighbour(double* image, int rows, int cols,
double r, double c, char mode,
double cval=*)
cdef inline double nearest_neighbour_interpolation(double* image, int rows,
int cols, double r,
double c, char mode,
double cval) nogil
cdef inline double bilinear_interpolation(double* image, int rows, int cols,
double r, double c, char mode,
double cval=*)
double cval) nogil
cdef inline double quadratic_interpolation(double x, double[3] f) nogil
cdef inline double biquadratic_interpolation(double* image, int rows, int cols,
double r, double c, char mode,
double cval) nogil
cdef inline double cubic_interpolation(double x, double[4] f) nogil
cdef inline double bicubic_interpolation(double* image, int rows, int cols,
double r, double c, char mode,
double cval) nogil
cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
char mode, double cval=*)
char mode, double cval) nogil
cdef inline int coord_map(int dim, int coord, char mode)
cdef inline int coord_map(int dim, int coord, char mode) nogil
skimage/_shared/interpolation.pyx
+189 -24
View File
@@ -5,24 +5,30 @@
from libc.math cimport ceil, floor, round
cdef inline double nearest_neighbour(double* image, int rows, int cols,
double r, double c, char mode,
double cval=0):
cdef inline double nearest_neighbour_interpolation(double* image, int rows,
int cols, double r,
double c, char mode,
double cval) nogil:
"""Nearest neighbour interpolation at a given position in the image.
Parameters
----------
image : double array
Input image.
rows, cols: int
rows, cols : int
Shape of image.
r, c : int
Position at which to interpolate.
mode : {'C', 'W', 'M'}
Wrapping mode. Constant, Wrap or Mirror.
mode : {'C', 'W', 'R', 'N'}
Wrapping mode. Constant, Wrap, Reflect or Nearest.
cval : double
Constant value to use for constant mode.
Returns
-------
value : double
Interpolated value.
"""
return get_pixel(image, rows, cols, <int>round(r), <int>round(c),
@@ -31,22 +37,27 @@ cdef inline double nearest_neighbour(double* image, int rows, int cols,
cdef inline double bilinear_interpolation(double* image, int rows, int cols,
double r, double c, char mode,
double cval=0):
double cval) nogil:
"""Bilinear interpolation at a given position in the image.
Parameters
----------
image : double array
Input image.
rows, cols: int
rows, cols : int
Shape of image.
r, c : int
Position at which to interpolate.
mode : {'C', 'W', 'M'}
Wrapping mode. Constant, Wrap or Mirror.
mode : {'C', 'W', 'R', 'N'}
Wrapping mode. Constant, Wrap, Reflect or Nearest.
cval : double
Constant value to use for constant mode.
Returns
-------
value : double
Interpolated value.
"""
cdef double dr, dc
cdef int minr, minc, maxr, maxc
@@ -64,23 +75,172 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols,
return (1 - dr) * top + dr * bottom
cdef inline double quadratic_interpolation(double x, double[3] f) nogil:
"""Quadratic interpolation.
Parameters
----------
x : double
Position in the interval [-1, 1].
f : double[4]
Function values at positions [-1, 0, 1].
Returns
-------
value : double
Interpolated value.
"""
return f[1] - 0.25 * (f[0] - f[2]) * x
cdef inline double biquadratic_interpolation(double* image, int rows, int cols,
double r, double c, char mode,
double cval) nogil:
"""Biquadratic interpolation at a given position in the image.
Parameters
----------
image : double array
Input image.
rows, cols : int
Shape of image.
r, c : int
Position at which to interpolate.
mode : {'C', 'W', 'R', 'N'}
Wrapping mode. Constant, Wrap, Reflect or Nearest.
cval : double
Constant value to use for constant mode.
Returns
-------
value : double
Interpolated value.
"""
cdef int r0 = <int>round(r)
cdef int c0 = <int>round(c)
if r < 0:
r0 -= 1
if c < 0:
c0 -= 1
# scale position to range [-1, 1]
cdef double xr = (r - r0) - 1
cdef double xc = (c - c0) - 1
if r == r0:
xr += 1
if c == c0:
xc += 1
cdef double fc[3], fr[3]
cdef int pr, pc
# row-wise cubic interpolation
for pr in range(r0, r0 + 3):
for pc in range(c0, c0 + 3):
fc[pc - c0] = get_pixel(image, rows, cols, pr, pc, mode, cval)
fr[pr - r0] = quadratic_interpolation(xc, fc)
# cubic interpolation for interpolated values of each row
return quadratic_interpolation(xr, fr)
cdef inline double cubic_interpolation(double x, double[4] f) nogil:
"""Cubic interpolation.
Parameters
----------
x : double
Position in the interval [0, 1].
f : double[4]
Function values at positions [0, 1/3, 2/3, 1].
Returns
-------
value : double
Interpolated value.
"""
return \
f[1] + 0.5 * x * \
(f[2] - f[0] + x * \
(2.0 * f[0] - 5.0 * f[1] + 4.0 * f[2] - f[3] + x * \
(3.0 * (f[1] - f[2]) + f[3] - f[0])))
cdef inline double bicubic_interpolation(double* image, int rows, int cols,
double r, double c, char mode,
double cval) nogil:
"""Bicubic interpolation at a given position in the image.
Parameters
----------
image : double array
Input image.
rows, cols : int
Shape of image.
r, c : int
Position at which to interpolate.
mode : {'C', 'W', 'R', 'N'}
Wrapping mode. Constant, Wrap, Reflect or Nearest.
cval : double
Constant value to use for constant mode.
Returns
-------
value : double
Interpolated value.
"""
cdef int r0 = <int>r - 1
cdef int c0 = <int>c - 1
if r < 0:
r0 -= 1
if c < 0:
c0 -= 1
# scale position to range [0, 1]
cdef double xr = (r - r0) / 3
cdef double xc = (c - c0) / 3
cdef double fc[4], fr[4]
cdef int pr, pc
# row-wise cubic interpolation
for pr in range(r0, r0 + 4):
for pc in range(c0, c0 + 4):
fc[pc - c0] = get_pixel(image, rows, cols, pr, pc, mode, cval)
fr[pr - r0] = cubic_interpolation(xc, fc)
# cubic interpolation for interpolated values of each row
return cubic_interpolation(xr, fr)
cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
char mode, double cval=0):
char mode, double cval) nogil:
"""Get a pixel from the image, taking wrapping mode into consideration.
Parameters
----------
image : double array
Input image.
rows, cols: int
rows, cols : int
Shape of image.
r, c : int
Position at which to get the pixel.
mode : {'C', 'W', 'M'}
Wrapping mode. Constant, Wrap or Mirror.
mode : {'C', 'W', 'R', 'N'}
Wrapping mode. Constant, Wrap, Reflect or Nearest.
cval : double
Constant value to use for constant mode.
Returns
-------
value : double
Pixel value at given position.
"""
if mode == 'C':
if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1):
@@ -91,7 +251,7 @@ cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
return image[coord_map(rows, r, mode) * cols + coord_map(cols, c, mode)]
cdef inline int coord_map(int dim, int coord, char mode):
cdef inline int coord_map(int dim, int coord, char mode) nogil:
"""
Wrap a coordinate, according to a given mode.
@@ -101,28 +261,33 @@ cdef inline int coord_map(int dim, int coord, char mode):
Maximum coordinate.
coord : int
Coord provided by user. May be < 0 or > dim.
mode : {'W', 'M'}
Whether to wrap or mirror the coordinate if it
mode : {'W', 'R', 'N'}
Whether to wrap or reflect the coordinate if it
falls outside [0, dim).
"""
dim = dim - 1
if mode == 'M': # mirror
if (coord < 0):
if mode == 'R': # reflect
if coord < 0:
# How many times times does the coordinate wrap?
if (<int>(-coord / dim) % 2 != 0):
if <int>(-coord / dim) % 2 != 0:
return dim - <int>(-coord % dim)
else:
return <int>(-coord % dim)
elif (coord > dim):
if (<int>(coord / dim) % 2 != 0):
elif coord > dim:
if <int>(coord / dim) % 2 != 0:
return <int>(dim - (coord % dim))
else:
return <int>(coord % dim)
elif mode == 'W': # wrap
if (coord < 0):
if coord < 0:
return <int>(dim - (-coord % dim))
elif (coord > dim):
elif coord > dim:
return <int>(coord % dim)
elif mode == 'N': # nearest
if coord < 0:
return 0
elif coord > dim:
return dim
return coord
skimage/feature/_texture.pyx
+1 -1
View File
@@ -132,7 +132,7 @@ def _local_binary_pattern(np.ndarray[double, ndim=2] image,
for c in range(image.shape[1]):
for i in range(P):
texture[i] = bilinear_interpolation(<double*>image.data,
rows, cols, r + coords[i, 0], c + coords[i, 1], 'C')
rows, cols, r + coords[i, 0], c + coords[i, 1], 'C', 0)
# signed / thresholded texture
for i in range(P):
if texture[i] - image[r, c] >= 0:
skimage/transform/__init__.py
+1 -2
View File
@@ -1,9 +1,8 @@
from .hough_transform import *
from .radon_transform import *
from .finite_radon_transform import *
from ._project import homography as fast_homography
from .integral import *
from ._geometric import (warp, warp_coords, estimate_transform,
SimilarityTransform, AffineTransform,
ProjectiveTransform, PolynomialTransform)
from ._warps import swirl, homography
from ._warps import resize, rotate, swirl, homography
skimage/transform/_geometric.py
+85 -59
View File
@@ -2,29 +2,7 @@ import math
import numpy as np
from scipy import ndimage
from skimage.util import img_as_float
def _stackcopy(a, b):
"""Copy b into each color layer of a, such that::
a[:,:,0] = a[:,:,1] = ... = b
Parameters
----------
a : (M, N) or (M, N, P) ndarray
Target array.
b : (M, N)
Source array.
Notes
-----
Color images are stored as an ``(M, N, 3)`` or ``(M, N, 4)`` arrays.
"""
if a.ndim == 3:
a[:] = b[:, :, np.newaxis]
else:
a[:] = b
from ._warps_cy import _warp_fast
class GeometricTransform(object):
@@ -603,12 +581,17 @@ class PolynomialTransform(GeometricTransform):
'then apply the forward transformation.')
TRANSFORMATIONS = {
TRANSFORMS = {
'similarity': SimilarityTransform,
'affine': AffineTransform,
'projective': ProjectiveTransform,
'polynomial': PolynomialTransform,
}
HOMOGRAPHY_TRANSFORMS = (
SimilarityTransform,
AffineTransform,
ProjectiveTransform
)
def estimate_transform(ttype, src, dst, **kwargs):
@@ -660,8 +643,8 @@ def estimate_transform(ttype, src, dst, **kwargs):
>>> warp(image, inverse_map=tform.inverse)
>>> # create transformation with explicit parameters
>>> tform2 = tf.SimilarityTransform()
>>> tform2.compose_implicit(scale=1.1, rotation=1, translation=(10, 20))
>>> tform2 = tf.SimilarityTransform(scale=1.1, rotation=1,
... translation=(10, 20))
>>> # unite transformations, applied in order from left to right
>>> tform3 = tform + tform2
@@ -669,11 +652,11 @@ def estimate_transform(ttype, src, dst, **kwargs):
"""
ttype = ttype.lower()
if ttype not in TRANSFORMATIONS:
if ttype not in TRANSFORMS:
raise ValueError('the transformation type \'%s\' is not'
'implemented' % ttype)
tform = TRANSFORMATIONS[ttype]()
tform = TRANSFORMS[ttype]()
tform.estimate(src, dst, **kwargs)
return tform
@@ -698,26 +681,44 @@ def matrix_transform(coords, matrix):
return ProjectiveTransform(matrix)(coords)
def warp_coords(orows, ocols, bands, coord_transform_fn,
dtype=np.float64):
def _stackcopy(a, b):
"""Copy b into each color layer of a, such that::
a[:,:,0] = a[:,:,1] = ... = b
Parameters
----------
a : (M, N) or (M, N, P) ndarray
Target array.
b : (M, N)
Source array.
Notes
-----
Color images are stored as an ``(M, N, 3)`` or ``(M, N, 4)`` arrays.
"""
if a.ndim == 3:
a[:] = b[:, :, np.newaxis]
else:
a[:] = b
def warp_coords(coord_map, shape, dtype=np.float64):
"""Build the source coordinates for the output pixels of an image warp.
Parameters
----------
orows : int
Number of output rows.
ocols : int
Number of output columns.
bands : int
Number of color bands (aka channels).
coord_transform_fn : callable like GeometricTransform.inverse
coord_map : callable like GeometricTransform.inverse
Return input coordinates for given output coordinates.
shape : tuple
Shape of output image ``(rows, cols[, bands])``.
dtype : np.dtype or string
dtype for return value (sane choices: float32 or float64)
dtype for return value (sane choices: float32 or float64).
Returns
-------
coords : (3, orows, ocols, bands) array of dtype `dtype`
coords : (ndim, rows, cols[, bands]) array of dtype `dtype`
Coordinates for `scipy.ndimage.map_coordinates`, that will yield
an image of shape (orows, ocols, bands) by drawing from source
points according to the `coord_transform_fn`.
@@ -745,23 +746,26 @@ def warp_coords(orows, ocols, bands, coord_transform_fn,
... xy[:, 0] -= 10
... return xy
>>>
>>> coords = warp_coords(30, 30, 3, shift_right)
>>> image = data.lena().astype(np.float32)
>>> coords = warp_coords(shift_right, image.shape)
>>> warped_image = map_coordinates(image, coords)
"""
coords = np.empty((3, orows, ocols, bands), dtype=dtype)
rows, cols = shape[0], shape[1]
coords_shape = [len(shape), rows, cols]
if len(shape) == 3:
coords_shape.append(shape[2])
coords = np.empty(coords_shape, dtype=dtype)
# Reshape grid coordinates into a (P, 2) array of (x, y) pairs
tf_coords = np.indices((ocols, orows), dtype=dtype).reshape(2, -1).T
tf_coords = np.indices((cols, rows), dtype=dtype).reshape(2, -1).T
# Map each (x, y) pair to the source image according to
# the user-provided mapping
tf_coords = coord_transform_fn(tf_coords)
tf_coords = coord_map(tf_coords)
# Reshape back to a (2, M, N) coordinate grid
tf_coords = tf_coords.T.reshape((-1, ocols, orows)).swapaxes(1, 2)
tf_coords = tf_coords.T.reshape((-1, cols, rows)).swapaxes(1, 2)
# Place the y-coordinate mapping
_stackcopy(coords[1, ...], tf_coords[0, ...])
@@ -769,8 +773,8 @@ def warp_coords(orows, ocols, bands, coord_transform_fn,
# Place the x-coordinate mapping
_stackcopy(coords[0, ...], tf_coords[1, ...])
# colour-coordinate mapping
coords[2, ...] = range(bands)
if len(shape) == 3:
coords[2, ...] = range(shape[2])
return coords
@@ -823,31 +827,53 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
if image.ndim < 2:
raise ValueError("Input must have more than 1 dimension.")
orig_ndim = image.ndim
image = np.atleast_3d(img_as_float(image))
ishape = np.array(image.shape)
bands = ishape[2]
if output_shape is None:
output_shape = ishape
# use fast Cython version for specific interpolation orders
if order in range(4) and not map_args:
matrix = None
if isinstance(inverse_map, HOMOGRAPHY_TRANSFORMS):
matrix = inverse_map._matrix
elif inverse_map.__name__ == 'inverse' \
and inverse_map.im_class in HOMOGRAPHY_TRANSFORMS:
matrix = np.linalg.inv(inverse_map.im_self._matrix)
if matrix is not None:
# transform all bands
dims = []
for dim in range(image.shape[2]):
dims.append(_warp_fast(image[..., dim], matrix,
output_shape=output_shape,
order=order, mode=mode, cval=cval))
out = np.dstack(dims)
if orig_ndim == 2:
out = out[..., 0]
rows, cols = output_shape[:2]
else: # use ndimage.map_coordinates
def coord_transform_fn(*args):
return inverse_map(*args, **map_args)
if output_shape is None:
output_shape = ishape
coords = warp_coords(rows, cols, bands, coord_transform_fn)
rows, cols = output_shape[:2]
# Prefilter not necessary for order 1 interpolation
prefilter = order > 1
mapped = ndimage.map_coordinates(image, coords, prefilter=prefilter,
mode=mode, order=order, cval=cval)
def coord_map(*args):
return inverse_map(*args, **map_args)
coords = warp_coords(coord_map, (rows, cols, bands))
# Prefilter not necessary for order 1 interpolation
prefilter = order > 1
out = ndimage.map_coordinates(image, coords, prefilter=prefilter,
mode=mode, order=order, cval=cval)
# The spline filters sometimes return results outside [0, 1],
# so clip to ensure valid data
clipped = np.clip(mapped, 0, 1)
clipped = np.clip(out, 0, 1)
if mode == 'constant' and not (0 <= cval <= 1):
clipped[mapped == cval] = cval
clipped[out == cval] = cval
# Remove singleton dim introduced by atleast_3d
return clipped.squeeze()
skimage/transform/_warps.py
+123 -2
View File
@@ -1,5 +1,119 @@
from ._geometric import warp, ProjectiveTransform
import numpy as np
from ._geometric import warp, SimilarityTransform, AffineTransform
def resize(image, output_shape, order=1, mode='constant', cval=0.):
"""Resize image.
Parameters
----------
image : ndarray
Input image.
output_shape : tuple or ndarray
Size of the generated output image `(rows, cols)`.
Returns
-------
resized : ndarray
Resized version of the input.
Other parameters
----------------
order : int
Order of splines used in interpolation. See
`scipy.ndimage.map_coordinates` for detail.
mode : string
How to handle values outside the image borders. See
`scipy.ndimage.map_coordinates` for detail.
cval : string
Used in conjunction with mode 'constant', the value outside
the image boundaries.
"""
rows, cols = output_shape
orig_rows, orig_cols = image.shape[0], image.shape[1]
rscale = float(orig_rows) / rows
cscale = float(orig_cols) / cols
# 3 control points necessary to estimate exact AffineTransform
src_corners = np.array([[1, 1], [1, rows], [cols, rows]]) - 1
dst_corners = np.zeros(src_corners.shape, dtype=np.double)
# take into account that 0th pixel is at position (0.5, 0.5)
dst_corners[:, 0] = cscale * (src_corners[:, 0] + 0.5) - 0.5
dst_corners[:, 1] = rscale * (src_corners[:, 1] + 0.5) - 0.5
tform = AffineTransform()
tform.estimate(src_corners, dst_corners)
return warp(image, tform, output_shape=output_shape, order=order,
mode=mode, cval=cval)
def rotate(image, angle, resize=False, order=1, mode='constant', cval=0.):
"""Rotate image by a certain angle around its center.
Parameters
----------
image : ndarray
Input image.
angle : float
Rotation angle in degrees in counter-clockwise direction.
resize: bool, optional
Determine whether the shape of the output image will be automatically
calculated, so the complete rotated image exactly fits. Default is
False.
Returns
-------
rotated : ndarray
Rotated version of the input.
Other parameters
----------------
order : int
Order of splines used in interpolation. See
`scipy.ndimage.map_coordinates` for detail.
mode : string
How to handle values outside the image borders. See
`scipy.ndimage.map_coordinates` for detail.
cval : string
Used in conjunction with mode 'constant', the value outside
the image boundaries.
"""
rows, cols = image.shape[0], image.shape[1]
# rotation around center
translation = np.array((cols, rows)) / 2. - 0.5
tform1 = SimilarityTransform(translation=-translation)
tform2 = SimilarityTransform(rotation=np.deg2rad(angle))
tform3 = SimilarityTransform(translation=translation)
tform = tform1 + tform2 + tform3
output_shape = None
if not resize:
# determine shape of output image
corners = np.array([[1, 1], [1, rows], [cols, rows], [cols, 1]])
corners = tform(corners - 1)
minc = corners[:, 0].min()
minr = corners[:, 1].min()
maxc = corners[:, 0].max()
maxr = corners[:, 1].max()
out_rows = maxr - minr + 1
out_cols = maxc - minc + 1
output_shape = np.ceil((out_rows, out_cols))
# fit output image in new shape
translation = ((cols - out_cols) / 2., (rows - out_rows) / 2.)
tform4 = SimilarityTransform(translation=translation)
tform = tform4 + tform
return warp(image, tform, output_shape=output_shape, order=order,
mode=mode, cval=cval)
def _swirl_mapping(xy, center, rotation, strength, radius):
x, y = xy.T
@@ -74,7 +188,14 @@ def swirl(image, center=None, strength=1, radius=100, rotation=0,
def homography(image, H, output_shape=None, order=1,
mode='constant', cval=0.):
"""Perform a projective transformation (homography) on an image.
"""
.. note:: Deprecated in skimage 0.7
`homography` will be removed in skimage 0.8, it is replaced by
`warp` because the latter provides the same functionality::
warp(image, ProjectiveTransform(H))
Perform a projective transformation (homography) on an image.
For each pixel, given its homogeneous coordinate :math:`\mathbf{x}
= [x, y, 1]^T`, its target position is calculated by multiplying
skimage/transform/_project.pyx → skimage/transform/_warps_cy.pyx
+34 -26
View File
@@ -5,12 +5,15 @@
cimport numpy as np
import numpy as np
from skimage._shared.interpolation cimport (nearest_neighbour,
bilinear_interpolation)
from cython.parallel import prange
from skimage._shared.interpolation cimport (nearest_neighbour_interpolation,
bilinear_interpolation,
biquadratic_interpolation,
bicubic_interpolation)
cdef inline _matrix_transform(double x, double y, double* H, double *x_,
double *y_):
cdef inline void _matrix_transform(double x, double y, double* H, double *x_,
double *y_) nogil:
"""Apply a homography to a coordinate.
Parameters
@@ -33,10 +36,9 @@ cdef inline _matrix_transform(double x, double y, double* H, double *x_,
y_[0] = yy / zz
def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
mode='constant', double cval=0):
"""
Projective transformation (homography).
"""Projective transformation (homography).
Perform a projective transformation (homography) of a
floating point image, using bi-linear interpolation.
@@ -72,7 +74,9 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
Order of interpolation::
* 0: Nearest-neighbour interpolation.
* 1: Bilinear interpolation (default).
mode : {'constant', 'mirror', 'wrap'}
* 2: Biquadratic interpolation (default).
* 3: Bicubic interpolation.
mode : {'constant', 'reflect', 'wrap'}
How to handle values outside the image borders.
cval : string
Used in conjunction with mode 'C' (constant), the value
@@ -83,18 +87,12 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] img = \
np.ascontiguousarray(image, dtype=np.double)
cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] M = \
np.ascontiguousarray(np.linalg.inv(H))
np.ascontiguousarray(H)
if mode not in ('constant', 'wrap', 'mirror'):
if mode not in ('constant', 'wrap', 'reflect', 'nearest'):
raise ValueError("Invalid mode specified. Please use "
"`constant`, `wrap` or `mirror`.")
cdef char mode_c
if mode == 'constant':
mode_c = ord('C')
elif mode == 'wrap':
mode_c = ord('W')
elif mode == 'mirror':
mode_c = ord('M')
"`constant`, `nearest`, `wrap` or `reflect`.")
cdef char mode_c = ord(mode[0].upper())
cdef int out_r, out_c
if output_shape is None:
@@ -104,22 +102,32 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1,
out_r = output_shape[0]
out_c = output_shape[1]
cdef np.ndarray[dtype=np.double_t, ndim=2] out = \
cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] out = \
np.zeros((out_r, out_c), dtype=np.double)
cdef double* out_data = <double*>out.data
cdef int tfr, tfc
cdef double r, c
cdef int rows = img.shape[0]
cdef int cols = img.shape[1]
for tfr in range(out_r):
cdef double (*interp_func)(double*, int, int, double, double,
char, double) nogil
if order == 0:
interp_func = nearest_neighbour_interpolation
elif order == 1:
interp_func = bilinear_interpolation
elif order == 2:
interp_func = biquadratic_interpolation
elif order == 3:
interp_func = bicubic_interpolation
for tfr in prange(out_r, nogil=True):
# make r, c thread local variables
r = c = 0
for tfc in range(out_c):
_matrix_transform(tfc, tfr, <double*>M.data, &c, &r)
if order == 0:
out[tfr, tfc] = nearest_neighbour(<double*>img.data, rows,
cols, r, c, mode_c)
elif order == 1:
out[tfr, tfc] = bilinear_interpolation(<double*>img.data, rows,
cols, r, c, mode_c)
out_data[tfr * out_r + tfc] = interp_func(<double*>img.data, rows,
cols, r, c, mode_c, cval)
return out
skimage/transform/radon_transform.py
+3 -3
View File
@@ -15,7 +15,7 @@ References:
from __future__ import division
import numpy as np
from scipy.fftpack import fftshift, fft, ifft
from ._project import homography
from ._warps_cy import _warp_fast
__all__ = ["radon", "iradon"]
@@ -77,8 +77,8 @@ def radon(image, theta=None):
return shift1.dot(R).dot(shift0)
for i in range(len(theta)):
rotated = homography(padded_image,
build_rotation(-theta[i]))
rotated = _warp_fast(padded_image,
np.linalg.inv(build_rotation(-theta[i])))
out[:, i] = rotated.sum(0)[::-1]
skimage/transform/setup.py
+5 -3
View File
@@ -14,13 +14,15 @@ def configuration(parent_package='', top_path=None):
config.add_data_dir('tests')
cython(['_hough_transform.pyx'], working_path=base_path)
cython(['_project.pyx'], working_path=base_path)
cython(['_warps_cy.pyx'], working_path=base_path)
config.add_extension('_hough_transform', sources=['_hough_transform.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('_project', sources=['_project.c'],
include_dirs=[get_numpy_include_dirs(), '../_shared'])
config.add_extension('_warps_cy', sources=['_warps_cy.c'],
include_dirs=[get_numpy_include_dirs(), '../_shared'],
extra_compile_args=['-fopenmp'],
extra_link_args=['-fopenmp'])
return config
skimage/transform/tests/test_geometric.py
+3 -4
View File
@@ -1,10 +1,9 @@
import numpy as np
from numpy.testing import assert_equal, assert_array_almost_equal
from skimage.transform._geometric import _stackcopy
from skimage.transform import (estimate_transform, SimilarityTransform,
AffineTransform, ProjectiveTransform,
PolynomialTransform)
from skimage.transform import (estimate_transform,
SimilarityTransform, AffineTransform,
ProjectiveTransform, PolynomialTransform)
SRC = np.array([
skimage/transform/tests/test_warps.py
+26 -8
View File
@@ -2,7 +2,7 @@ from numpy.testing import assert_array_almost_equal, run_module_suite
import numpy as np
from scipy.ndimage import map_coordinates
from skimage.transform import (warp, warp_coords, fast_homography,
from skimage.transform import (warp, warp_coords, rotate, resize,
AffineTransform,
ProjectiveTransform,
SimilarityTransform)
@@ -55,11 +55,12 @@ def test_fast_homography():
H[:2, 2] = [tx, ty]
tform = ProjectiveTransform(H)
coords = warp_coords(tform.inverse, (img.shape[0], img.shape[1]))
for order in range(2):
for mode in ('constant', 'mirror', 'wrap'):
p0 = warp(img, tform.inverse, mode=mode, order=order)
p1 = fast_homography(img, H, mode=mode, order=order)
for order in range(4):
for mode in ('constant', 'reflect', 'wrap', 'nearest'):
p0 = map_coordinates(img, coords, mode=mode, order=order)
p1 = warp(img, tform, mode=mode, order=order)
# import matplotlib.pyplot as plt
# f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
@@ -73,6 +74,22 @@ def test_fast_homography():
assert d < 0.001
def test_rotate():
x = np.zeros((5, 5), dtype=np.double)
x[1, 1] = 1
x90 = rotate(x, 90)
assert_array_almost_equal(x90, np.rot90(x))
def test_resize():
x = np.zeros((5, 5), dtype=np.double)
x[1, 1] = 1
resized = resize(x, (10, 10), order=0)
ref = np.zeros((10, 10))
ref[2:4, 2:4] = 1
assert_array_almost_equal(resized, ref)
def test_swirl():
image = img_as_float(data.checkerboard())
@@ -85,8 +102,9 @@ def test_swirl():
def test_const_cval_out_of_range():
img = np.random.randn(100, 100)
warped = warp(img, AffineTransform(translation=(10, 10)), cval=-10)
assert np.sum(warped < 0) == (2 * 100 * 10 - 10 * 10)
cval = - 10
warped = warp(img, AffineTransform(translation=(10, 10)), cval=cval)
assert np.sum(warped == cval) == (2 * 100 * 10 - 10 * 10)
def test_warp_identity():
@@ -107,7 +125,7 @@ def test_warp_coords_example():
image = data.lena().astype(np.float32)
assert 3 == image.shape[2]
tform = SimilarityTransform(translation=(0, -10))
coords = warp_coords(30, 30, 3, tform)
coords = warp_coords(tform, (30, 30, 3))
warped_image1 = map_coordinates(image[:, :, 0], coords[:2])