From 77f1e0ba473240fb9fcccd5d9ec884b9556b63c4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 26 Aug 2012 08:52:20 +0200 Subject: [PATCH 01/24] Add deprecation warning in doc string to homography --- skimage/transform/_project.pyx | 3 +-- skimage/transform/_warps.py | 6 +++++- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/skimage/transform/_project.pyx b/skimage/transform/_project.pyx index 4df6e001..201785ad 100644 --- a/skimage/transform/_project.pyx +++ b/skimage/transform/_project.pyx @@ -35,8 +35,7 @@ cdef inline _matrix_transform(double x, double y, double* H, double *x_, def homography(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. diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py index f09f7944..4c499d6b 100644 --- a/skimage/transform/_warps.py +++ b/skimage/transform/_warps.py @@ -74,7 +74,11 @@ 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. + """ + .. deprecated:: + 0.7 + + 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 From a6532a8dae3f7e4a211a8f5825a4110f0cf5dfd0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 26 Aug 2012 15:39:48 +0200 Subject: [PATCH 02/24] Refactor image warps * Fix cval bug in interpolation which was ignored * Remove fast_homography as standalone function and automatically include functionality in warp * Fix bug in warp_coords for graylevel images * move warp functions to warp file --- skimage/_shared/interpolation.pyx | 18 +- skimage/transform/__init__.py | 5 +- skimage/transform/_geometric.py | 186 +-------------- skimage/transform/_warps.py | 212 +++++++++++++++++- .../transform/{_project.pyx => _warps_cy.pyx} | 18 +- skimage/transform/setup.py | 4 +- skimage/transform/tests/test_warps.py | 16 +- 7 files changed, 245 insertions(+), 214 deletions(-) rename skimage/transform/{_project.pyx => _warps_cy.pyx} (90%) diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx index 137004e5..83722aba 100644 --- a/skimage/_shared/interpolation.pyx +++ b/skimage/_shared/interpolation.pyx @@ -18,8 +18,8 @@ cdef inline double nearest_neighbour(double* image, int rows, int cols, 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'} + Wrapping mode. Constant, Wrap or Reflect. cval : double Constant value to use for constant mode. @@ -42,8 +42,8 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols, 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'} + Wrapping mode. Constant, Wrap or Reflect. cval : double Constant value to use for constant mode. @@ -76,8 +76,8 @@ cdef inline double get_pixel(double* image, int rows, int cols, int r, int c, 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'} + Wrapping mode. Constant, Wrap or Reflect. cval : double Constant value to use for constant mode. @@ -101,13 +101,13 @@ 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'} + Whether to wrap or reflect the coordinate if it falls outside [0, dim). """ dim = dim - 1 - if mode == 'M': # mirror + if mode == 'R': # reflect if (coord < 0): # How many times times does the coordinate wrap? if ((-coord / dim) % 2 != 0): diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index b0eee512..0735267a 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -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, +from ._geometric import (estimate_transform, SimilarityTransform, AffineTransform, ProjectiveTransform, PolynomialTransform) -from ._warps import swirl, homography +from ._warps import warp, warp_coords, swirl, homography diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index c3b47a3b..e4799375 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -1,30 +1,5 @@ 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 class GeometricTransform(object): @@ -603,7 +578,7 @@ class PolynomialTransform(GeometricTransform): 'then apply the forward transformation.') -TRANSFORMATIONS = { +TRANSFORMS = { 'similarity': SimilarityTransform, 'affine': AffineTransform, 'projective': ProjectiveTransform, @@ -669,11 +644,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 @@ -696,158 +671,3 @@ def matrix_transform(coords, matrix): """ return ProjectiveTransform(matrix)(coords) - - -def warp_coords(orows, ocols, bands, coord_transform_fn, - 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 - Return input coordinates for given output coordinates. - dtype : np.dtype or string - dtype for return value (sane choices: float32 or float64) - - Returns - ------- - coords : (3, orows, ocols, 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`. - - Notes - ----- - This is a lower-level routine that produces the source coordinates used by - `warp()`. - - It is provided separately from `warp` to give additional flexibility to - users who would like, for example, to re-use a particular coordinate - mapping, to use specific dtypes at various points along the the - image-warping process, or to implement different post-processing logic - than `warp` performs after the call to `ndimage.map_coordinates`. - - - Examples - -------- - Produce a coordinate map that Shifts an image to the right: - - >>> from skimage import data - >>> from scipy.ndimage import map_coordinates - >>> - >>> def shift_right(xy): - ... xy[:, 0] -= 10 - ... return xy - >>> - >>> coords = warp_coords(30, 30, 3, shift_right) - >>> image = data.lena().astype(np.float32) - >>> warped_image = map_coordinates(image, coords) - - """ - - coords = np.empty((3, orows, ocols, bands), 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 - - # Map each (x, y) pair to the source image according to - # the user-provided mapping - tf_coords = coord_transform_fn(tf_coords) - - # Reshape back to a (2, M, N) coordinate grid - tf_coords = tf_coords.T.reshape((-1, ocols, orows)).swapaxes(1, 2) - - # Place the y-coordinate mapping - _stackcopy(coords[1, ...], tf_coords[0, ...]) - - # Place the x-coordinate mapping - _stackcopy(coords[0, ...], tf_coords[1, ...]) - - # colour-coordinate mapping - coords[2, ...] = range(bands) - - return coords - - -def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, - mode='constant', cval=0., reverse_map=None): - """Warp an image according to a given coordinate transformation. - - Parameters - ---------- - image : 2-D array - Input image. - inverse_map : transformation object, callable xy = f(xy, **kwargs) - Inverse coordinate map. A function that transforms a (N, 2) array of - ``(x, y)`` coordinates in the *output image* into their corresponding - coordinates in the *source image* (e.g. a transformation object or its - inverse). - map_args : dict, optional - Keyword arguments passed to `inverse_map`. - output_shape : tuple (rows, cols) - Shape of the output image generated. - 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 : float - Used in conjunction with mode 'constant', the value outside - the image boundaries. - - Examples - -------- - Shift an image to the right: - - >>> from skimage import data - >>> image = data.camera() - >>> - >>> def shift_right(xy): - ... xy[:, 0] -= 10 - ... return xy - >>> - >>> warp(image, shift_right) - - """ - # Backward API compatibility - if reverse_map is not None: - inverse_map = reverse_map - - if image.ndim < 2: - raise ValueError("Input must have more than 1 dimension.") - - image = np.atleast_3d(img_as_float(image)) - ishape = np.array(image.shape) - bands = ishape[2] - - if output_shape is None: - output_shape = ishape - - rows, cols = output_shape[:2] - - def coord_transform_fn(*args): - return inverse_map(*args, **map_args) - - coords = warp_coords(rows, cols, bands, coord_transform_fn) - - # Prefilter not necessary for order 1 interpolation - prefilter = order > 1 - mapped = 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) - - if mode == 'constant' and not (0 <= cval <= 1): - clipped[mapped == cval] = cval - - # Remove singleton dim introduced by atleast_3d - return clipped.squeeze() diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py index 4c499d6b..db7e5eba 100644 --- a/skimage/transform/_warps.py +++ b/skimage/transform/_warps.py @@ -1,5 +1,215 @@ -from ._geometric import warp, ProjectiveTransform import numpy as np +from scipy import ndimage +from skimage.util import img_as_float +from ._geometric import (SimilarityTransform, AffineTransform, + ProjectiveTransform) +from ._warps_cy import _warp_fast + + +HOMOGRAPHY_TRANSFORMS = ( + SimilarityTransform, + AffineTransform, + ProjectiveTransform +) + + +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 + ---------- + 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). + + Returns + ------- + 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`. + + Notes + ----- + This is a lower-level routine that produces the source coordinates used by + `warp()`. + + It is provided separately from `warp` to give additional flexibility to + users who would like, for example, to re-use a particular coordinate + mapping, to use specific dtypes at various points along the the + image-warping process, or to implement different post-processing logic + than `warp` performs after the call to `ndimage.map_coordinates`. + + + Examples + -------- + Produce a coordinate map that Shifts an image to the right: + + >>> from skimage import data + >>> from scipy.ndimage import map_coordinates + >>> + >>> def shift_right(xy): + ... xy[:, 0] -= 10 + ... return xy + >>> + >>> coords = warp_coords(30, 30, 3, shift_right) + >>> image = data.lena().astype(np.float32) + >>> warped_image = map_coordinates(image, coords) + + """ + 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((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_map(tf_coords) + + # Reshape back to a (2, M, N) coordinate grid + tf_coords = tf_coords.T.reshape((-1, cols, rows)).swapaxes(1, 2) + + # Place the y-coordinate mapping + _stackcopy(coords[1, ...], tf_coords[0, ...]) + + # Place the x-coordinate mapping + _stackcopy(coords[0, ...], tf_coords[1, ...]) + + if len(shape) == 3: + coords[2, ...] = range(shape[2]) + + return coords + + +def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, + mode='constant', cval=0., reverse_map=None): + """Warp an image according to a given coordinate transformation. + + Parameters + ---------- + image : 2-D array + Input image. + inverse_map : transformation object, callable xy = f(xy, **kwargs) + Inverse coordinate map. A function that transforms a (N, 2) array of + ``(x, y)`` coordinates in the *output image* into their corresponding + coordinates in the *source image* (e.g. a transformation object or its + inverse). + map_args : dict, optional + Keyword arguments passed to `inverse_map`. + output_shape : tuple (rows, cols) + Shape of the output image generated. + 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 : float + Used in conjunction with mode 'constant', the value outside + the image boundaries. + + Examples + -------- + Shift an image to the right: + + >>> from skimage import data + >>> image = data.camera() + >>> + >>> def shift_right(xy): + ... xy[:, 0] -= 10 + ... return xy + >>> + >>> warp(image, shift_right) + + """ + # Backward API compatibility + if reverse_map is not None: + inverse_map = reverse_map + + 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] + + # use fast Cython version for specific parameters + fast_modes = ('constant', 'reflect', 'wrap') + if order in (0, 1) and mode in fast_modes 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] + return out + + if output_shape is None: + output_shape = ishape + + rows, cols = output_shape[:2] + + 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 + mapped = 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) + + if mode == 'constant' and not (0 <= cval <= 1): + clipped[mapped == cval] = cval + + # Remove singleton dim introduced by atleast_3d + return clipped.squeeze() + def _swirl_mapping(xy, center, rotation, strength, radius): x, y = xy.T diff --git a/skimage/transform/_project.pyx b/skimage/transform/_warps_cy.pyx similarity index 90% rename from skimage/transform/_project.pyx rename to skimage/transform/_warps_cy.pyx index 201785ad..68ccef27 100644 --- a/skimage/transform/_project.pyx +++ b/skimage/transform/_warps_cy.pyx @@ -33,7 +33,7 @@ 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). @@ -71,7 +71,7 @@ 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'} + mode : {'constant', 'reflect', 'wrap'} How to handle values outside the image borders. cval : string Used in conjunction with mode 'C' (constant), the value @@ -82,18 +82,18 @@ 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'): raise ValueError("Invalid mode specified. Please use " - "`constant`, `wrap` or `mirror`.") + "`constant`, `wrap` or `reflect`.") 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') + elif mode == 'reflect': + mode_c = ord('R') cdef int out_r, out_c if output_shape is None: @@ -116,9 +116,9 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None, int order=1, _matrix_transform(tfc, tfr, M.data, &c, &r) if order == 0: out[tfr, tfc] = nearest_neighbour(img.data, rows, - cols, r, c, mode_c) + cols, r, c, mode_c, cval) elif order == 1: out[tfr, tfc] = bilinear_interpolation(img.data, rows, - cols, r, c, mode_c) + cols, r, c, mode_c, cval) return out diff --git a/skimage/transform/setup.py b/skimage/transform/setup.py index 75210b54..0e415bad 100644 --- a/skimage/transform/setup.py +++ b/skimage/transform/setup.py @@ -14,12 +14,12 @@ 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'], + config.add_extension('_warps_cy', sources=['_warps_cy.c'], include_dirs=[get_numpy_include_dirs(), '../_shared']) return config diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py index 7be3151b..f896a452 100644 --- a/skimage/transform/tests/test_warps.py +++ b/skimage/transform/tests/test_warps.py @@ -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, 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 mode in ('constant', 'reflect', 'wrap'): + 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) @@ -85,8 +86,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 +109,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]) From cfea01b9fffba5f1b217bdae92865102373e6d65 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 26 Aug 2012 15:51:09 +0200 Subject: [PATCH 03/24] Fix import error in geometric test cases --- skimage/transform/tests/test_geometric.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index abe8b2e7..f035631d 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -1,7 +1,7 @@ import numpy as np from numpy.testing import assert_equal, assert_array_almost_equal -from skimage.transform._geometric import _stackcopy +from skimage.transform._warps import _stackcopy from skimage.transform import (estimate_transform, SimilarityTransform, AffineTransform, ProjectiveTransform, PolynomialTransform) From 772a1cb4b0929093119ca166ac79b4eeb2856392 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 26 Aug 2012 17:31:22 +0200 Subject: [PATCH 04/24] Update radon transform with new warp function --- skimage/transform/radon_transform.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py index f80cc3cb..0213b2bc 100644 --- a/skimage/transform/radon_transform.py +++ b/skimage/transform/radon_transform.py @@ -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] From a0649791aed1f2a75657cd94c29fdd4f5119c57d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 26 Aug 2012 17:33:23 +0200 Subject: [PATCH 05/24] Fix doc string example of warp_coords --- skimage/transform/_warps.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py index db7e5eba..e796b454 100644 --- a/skimage/transform/_warps.py +++ b/skimage/transform/_warps.py @@ -78,8 +78,8 @@ def warp_coords(coord_map, shape, dtype=np.float64): ... 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) """ From 6f4c2ad2681e4c9b917a26dbe4305a4ef5ceedf6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 26 Aug 2012 18:46:13 +0200 Subject: [PATCH 06/24] Add function for image rotation --- skimage/transform/_warps.py | 62 +++++++++++++++++++++++++++ skimage/transform/tests/test_warps.py | 9 +++- 2 files changed, 70 insertions(+), 1 deletion(-) diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py index e796b454..7b0fbeed 100644 --- a/skimage/transform/_warps.py +++ b/skimage/transform/_warps.py @@ -211,6 +211,68 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, return clipped.squeeze() +def rotate(image, angle, preserve_shape=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. + preserve_shape : 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. + tform1 = SimilarityTransform(translation=-translation) + tform2 = SimilarityTransform(rotation=np.deg2rad(angle)) + tform3 = SimilarityTransform(translation=translation) + tform = tform1 + tform2 + tform3 + + output_shape = None + if not preserve_shape: + # determine shape of output image + corners = tform([[0, 0], [0, rows], [cols, 0], [cols, rows]]) + corners = np.round(corners, 4) + minc = np.floor(corners[:, 0].min()) + minr = np.floor(corners[:, 1].min()) + maxc = np.ceil(corners[:, 0].max()) + maxr = np.ceil(corners[:, 1].max()) + output_shape = [maxr - minr, maxc - minc] + + # fit output image in new shape + tform4 = SimilarityTransform(translation=(minc, minr + 1)) + 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 x0, y0 = center diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py index f896a452..d5ffff59 100644 --- a/skimage/transform/tests/test_warps.py +++ b/skimage/transform/tests/test_warps.py @@ -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, +from skimage.transform import (warp, warp_coords, rotate, AffineTransform, ProjectiveTransform, SimilarityTransform) @@ -74,6 +74,13 @@ 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_swirl(): image = img_as_float(data.checkerboard()) From 6ec4c21cf7af09401aabadff79898fe783efe9bd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 26 Aug 2012 22:55:35 +0200 Subject: [PATCH 07/24] Fix import of rotate function --- skimage/transform/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index 0735267a..1ee519f2 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -5,4 +5,4 @@ from .integral import * from ._geometric import (estimate_transform, SimilarityTransform, AffineTransform, ProjectiveTransform, PolynomialTransform) -from ._warps import warp, warp_coords, swirl, homography +from ._warps import warp, warp_coords, rotate, swirl, homography From 00e1d14294ff9c11e06f49e8d366b4b60ad35493 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 26 Aug 2012 22:56:57 +0200 Subject: [PATCH 08/24] Fix rotate translation bug --- skimage/transform/_warps.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py index 7b0fbeed..e9cb61e5 100644 --- a/skimage/transform/_warps.py +++ b/skimage/transform/_warps.py @@ -248,7 +248,7 @@ def rotate(image, angle, preserve_shape=False, order=1, rows, cols = image.shape[0], image.shape[1] # rotation around center - translation = np.array((cols, rows)) / 2. + translation = np.floor(np.array((cols, rows )) / 2.) tform1 = SimilarityTransform(translation=-translation) tform2 = SimilarityTransform(rotation=np.deg2rad(angle)) tform3 = SimilarityTransform(translation=translation) @@ -257,16 +257,19 @@ def rotate(image, angle, preserve_shape=False, order=1, output_shape = None if not preserve_shape: # determine shape of output image - corners = tform([[0, 0], [0, rows], [cols, 0], [cols, rows]]) - corners = np.round(corners, 4) - minc = np.floor(corners[:, 0].min()) - minr = np.floor(corners[:, 1].min()) - maxc = np.ceil(corners[:, 0].max()) - maxr = np.ceil(corners[:, 1].max()) - output_shape = [maxr - minr, maxc - minc] + corners = np.array([[1, 1], [1, rows], [cols, 1], [cols, rows]]) + corners = tform2(tform1(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 = (out_rows, out_cols) # fit output image in new shape - tform4 = SimilarityTransform(translation=(minc, minr + 1)) + 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, From c0538c03e90cda78670855590618eb113ab23824 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 27 Aug 2012 13:21:11 +0200 Subject: [PATCH 09/24] Fix translation bug in image rotate --- skimage/transform/_warps.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py index e9cb61e5..14f0e187 100644 --- a/skimage/transform/_warps.py +++ b/skimage/transform/_warps.py @@ -248,7 +248,7 @@ def rotate(image, angle, preserve_shape=False, order=1, rows, cols = image.shape[0], image.shape[1] # rotation around center - translation = np.floor(np.array((cols, rows )) / 2.) + translation = np.array((cols, rows)) / 2. - 0.5 tform1 = SimilarityTransform(translation=-translation) tform2 = SimilarityTransform(rotation=np.deg2rad(angle)) tform3 = SimilarityTransform(translation=translation) @@ -257,15 +257,15 @@ def rotate(image, angle, preserve_shape=False, order=1, output_shape = None if not preserve_shape: # determine shape of output image - corners = np.array([[1, 1], [1, rows], [cols, 1], [cols, rows]]) - corners = tform2(tform1(corners - 1)) + 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 = (out_rows, out_cols) + output_shape = np.ceil((out_rows, out_cols)) # fit output image in new shape translation = ((cols - out_cols) / 2., (rows - out_rows) / 2.) From f4a3b44355e46be1c69feae9a57a1e0899ef8244 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 27 Aug 2012 13:23:43 +0200 Subject: [PATCH 10/24] Rename parameter of image rotation --- skimage/transform/_warps.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py index 14f0e187..2a2d7828 100644 --- a/skimage/transform/_warps.py +++ b/skimage/transform/_warps.py @@ -211,8 +211,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, return clipped.squeeze() -def rotate(image, angle, preserve_shape=False, order=1, - mode='constant', cval=0.): +def rotate(image, angle, resize=False, order=1, mode='constant', cval=0.): """Rotate image by a certain angle around its center. Parameters @@ -221,7 +220,7 @@ def rotate(image, angle, preserve_shape=False, order=1, Input image. angle : float Rotation angle in degrees in counter-clockwise direction. - preserve_shape : bool, optional + 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. @@ -255,7 +254,7 @@ def rotate(image, angle, preserve_shape=False, order=1, tform = tform1 + tform2 + tform3 output_shape = None - if not preserve_shape: + if not resize: # determine shape of output image corners = np.array([[1, 1], [1, rows], [cols, rows], [cols, 1]]) corners = tform(corners - 1) From 9a46111e761e527d337a628f1c263dfebe89cb2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 27 Aug 2012 13:35:36 +0200 Subject: [PATCH 11/24] Apply numpy doc style for deprecation warning --- skimage/transform/_warps.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py index 2a2d7828..c6e22a19 100644 --- a/skimage/transform/_warps.py +++ b/skimage/transform/_warps.py @@ -349,8 +349,11 @@ def swirl(image, center=None, strength=1, radius=100, rotation=0, def homography(image, H, output_shape=None, order=1, mode='constant', cval=0.): """ - .. deprecated:: - 0.7 + .. 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. From 501c401ddad3b2614abd99d3513a8c19fea25eca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 27 Aug 2012 18:05:18 +0200 Subject: [PATCH 12/24] Move image transform functions to _geometric file --- skimage/transform/__init__.py | 4 +- skimage/transform/_geometric.py | 206 +++++++++++++++++++++ skimage/transform/_warps.py | 211 +--------------------- skimage/transform/tests/test_geometric.py | 9 +- 4 files changed, 213 insertions(+), 217 deletions(-) diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index 1ee519f2..511c8738 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -2,7 +2,7 @@ from .hough_transform import * from .radon_transform import * from .finite_radon_transform import * from .integral import * -from ._geometric import (estimate_transform, +from ._geometric import (warp, warp_coords, estimate_transform, SimilarityTransform, AffineTransform, ProjectiveTransform, PolynomialTransform) -from ._warps import warp, warp_coords, rotate, swirl, homography +from ._warps import rotate, swirl, homography diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index e4799375..f11f1b44 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -1,5 +1,8 @@ import math import numpy as np +from scipy import ndimage +from skimage.util import img_as_float +from ._warps_cy import _warp_fast class GeometricTransform(object): @@ -584,6 +587,11 @@ TRANSFORMS = { 'projective': ProjectiveTransform, 'polynomial': PolynomialTransform, } +HOMOGRAPHY_TRANSFORMS = ( + SimilarityTransform, + AffineTransform, + ProjectiveTransform +) def estimate_transform(ttype, src, dst, **kwargs): @@ -671,3 +679,201 @@ def matrix_transform(coords, matrix): """ return ProjectiveTransform(matrix)(coords) + + +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 + ---------- + 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). + + Returns + ------- + 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`. + + Notes + ----- + This is a lower-level routine that produces the source coordinates used by + `warp()`. + + It is provided separately from `warp` to give additional flexibility to + users who would like, for example, to re-use a particular coordinate + mapping, to use specific dtypes at various points along the the + image-warping process, or to implement different post-processing logic + than `warp` performs after the call to `ndimage.map_coordinates`. + + + Examples + -------- + Produce a coordinate map that Shifts an image to the right: + + >>> from skimage import data + >>> from scipy.ndimage import map_coordinates + >>> + >>> def shift_right(xy): + ... xy[:, 0] -= 10 + ... return xy + >>> + >>> image = data.lena().astype(np.float32) + >>> coords = warp_coords(shift_right, image.shape) + >>> warped_image = map_coordinates(image, coords) + + """ + 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((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_map(tf_coords) + + # Reshape back to a (2, M, N) coordinate grid + tf_coords = tf_coords.T.reshape((-1, cols, rows)).swapaxes(1, 2) + + # Place the y-coordinate mapping + _stackcopy(coords[1, ...], tf_coords[0, ...]) + + # Place the x-coordinate mapping + _stackcopy(coords[0, ...], tf_coords[1, ...]) + + if len(shape) == 3: + coords[2, ...] = range(shape[2]) + + return coords + + +def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, + mode='constant', cval=0., reverse_map=None): + """Warp an image according to a given coordinate transformation. + + Parameters + ---------- + image : 2-D array + Input image. + inverse_map : transformation object, callable xy = f(xy, **kwargs) + Inverse coordinate map. A function that transforms a (N, 2) array of + ``(x, y)`` coordinates in the *output image* into their corresponding + coordinates in the *source image* (e.g. a transformation object or its + inverse). + map_args : dict, optional + Keyword arguments passed to `inverse_map`. + output_shape : tuple (rows, cols) + Shape of the output image generated. + 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 : float + Used in conjunction with mode 'constant', the value outside + the image boundaries. + + Examples + -------- + Shift an image to the right: + + >>> from skimage import data + >>> image = data.camera() + >>> + >>> def shift_right(xy): + ... xy[:, 0] -= 10 + ... return xy + >>> + >>> warp(image, shift_right) + + """ + # Backward API compatibility + if reverse_map is not None: + inverse_map = reverse_map + + 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] + + # use fast Cython version for specific parameters + fast_modes = ('constant', 'reflect', 'wrap') + if order in (0, 1) and mode in fast_modes 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] + return out + + if output_shape is None: + output_shape = ishape + + rows, cols = output_shape[:2] + + 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 + mapped = 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) + + if mode == 'constant' and not (0 <= cval <= 1): + clipped[mapped == cval] = cval + + # Remove singleton dim introduced by atleast_3d + return clipped.squeeze() diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py index c6e22a19..4e5a7245 100644 --- a/skimage/transform/_warps.py +++ b/skimage/transform/_warps.py @@ -1,214 +1,5 @@ import numpy as np -from scipy import ndimage -from skimage.util import img_as_float -from ._geometric import (SimilarityTransform, AffineTransform, - ProjectiveTransform) -from ._warps_cy import _warp_fast - - -HOMOGRAPHY_TRANSFORMS = ( - SimilarityTransform, - AffineTransform, - ProjectiveTransform -) - - -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 - ---------- - 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). - - Returns - ------- - 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`. - - Notes - ----- - This is a lower-level routine that produces the source coordinates used by - `warp()`. - - It is provided separately from `warp` to give additional flexibility to - users who would like, for example, to re-use a particular coordinate - mapping, to use specific dtypes at various points along the the - image-warping process, or to implement different post-processing logic - than `warp` performs after the call to `ndimage.map_coordinates`. - - - Examples - -------- - Produce a coordinate map that Shifts an image to the right: - - >>> from skimage import data - >>> from scipy.ndimage import map_coordinates - >>> - >>> def shift_right(xy): - ... xy[:, 0] -= 10 - ... return xy - >>> - >>> image = data.lena().astype(np.float32) - >>> coords = warp_coords(shift_right, image.shape) - >>> warped_image = map_coordinates(image, coords) - - """ - 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((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_map(tf_coords) - - # Reshape back to a (2, M, N) coordinate grid - tf_coords = tf_coords.T.reshape((-1, cols, rows)).swapaxes(1, 2) - - # Place the y-coordinate mapping - _stackcopy(coords[1, ...], tf_coords[0, ...]) - - # Place the x-coordinate mapping - _stackcopy(coords[0, ...], tf_coords[1, ...]) - - if len(shape) == 3: - coords[2, ...] = range(shape[2]) - - return coords - - -def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, - mode='constant', cval=0., reverse_map=None): - """Warp an image according to a given coordinate transformation. - - Parameters - ---------- - image : 2-D array - Input image. - inverse_map : transformation object, callable xy = f(xy, **kwargs) - Inverse coordinate map. A function that transforms a (N, 2) array of - ``(x, y)`` coordinates in the *output image* into their corresponding - coordinates in the *source image* (e.g. a transformation object or its - inverse). - map_args : dict, optional - Keyword arguments passed to `inverse_map`. - output_shape : tuple (rows, cols) - Shape of the output image generated. - 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 : float - Used in conjunction with mode 'constant', the value outside - the image boundaries. - - Examples - -------- - Shift an image to the right: - - >>> from skimage import data - >>> image = data.camera() - >>> - >>> def shift_right(xy): - ... xy[:, 0] -= 10 - ... return xy - >>> - >>> warp(image, shift_right) - - """ - # Backward API compatibility - if reverse_map is not None: - inverse_map = reverse_map - - 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] - - # use fast Cython version for specific parameters - fast_modes = ('constant', 'reflect', 'wrap') - if order in (0, 1) and mode in fast_modes 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] - return out - - if output_shape is None: - output_shape = ishape - - rows, cols = output_shape[:2] - - 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 - mapped = 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) - - if mode == 'constant' and not (0 <= cval <= 1): - clipped[mapped == cval] = cval - - # Remove singleton dim introduced by atleast_3d - return clipped.squeeze() +from ._geometric import warp, SimilarityTransform def rotate(image, angle, resize=False, order=1, mode='constant', cval=0.): diff --git a/skimage/transform/tests/test_geometric.py b/skimage/transform/tests/test_geometric.py index f035631d..57e95d46 100644 --- a/skimage/transform/tests/test_geometric.py +++ b/skimage/transform/tests/test_geometric.py @@ -1,10 +1,9 @@ import numpy as np from numpy.testing import assert_equal, assert_array_almost_equal - -from skimage.transform._warps import _stackcopy -from skimage.transform import (estimate_transform, SimilarityTransform, - AffineTransform, ProjectiveTransform, - PolynomialTransform) +from skimage.transform._geometric import _stackcopy +from skimage.transform import (estimate_transform, + SimilarityTransform, AffineTransform, + ProjectiveTransform, PolynomialTransform) SRC = np.array([ From 99e4264e15cd34d3875bd311470ce935a16c145b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Mon, 27 Aug 2012 18:15:39 +0200 Subject: [PATCH 13/24] Add missing return type to matrix transform function --- skimage/transform/_warps_cy.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx index 68ccef27..2851a728 100644 --- a/skimage/transform/_warps_cy.pyx +++ b/skimage/transform/_warps_cy.pyx @@ -9,8 +9,8 @@ from skimage._shared.interpolation cimport (nearest_neighbour, bilinear_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_): """Apply a homography to a coordinate. Parameters From 9b44e24f8ed9d711b8127ad8097acb9a0855304d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 30 Aug 2012 07:33:39 +0200 Subject: [PATCH 14/24] Use function pointer for different interpolation methods --- skimage/_shared/interpolation.pxd | 6 +++--- skimage/_shared/interpolation.pyx | 6 +++--- skimage/feature/_texture.pyx | 2 +- skimage/transform/_warps_cy.pyx | 15 +++++++++------ 4 files changed, 16 insertions(+), 13 deletions(-) diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd index c883d00d..3e8e74e9 100644 --- a/skimage/_shared/interpolation.pxd +++ b/skimage/_shared/interpolation.pxd @@ -1,13 +1,13 @@ cdef inline double nearest_neighbour(double* image, int rows, int cols, double r, double c, char mode, - double cval=*) + double cval) cdef inline double bilinear_interpolation(double* image, int rows, int cols, double r, double c, char mode, - double cval=*) + double cval) cdef inline double get_pixel(double* image, int rows, int cols, int r, int c, - char mode, double cval=*) + char mode, double cval) cdef inline int coord_map(int dim, int coord, char mode) diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx index 83722aba..63f5b70c 100644 --- a/skimage/_shared/interpolation.pyx +++ b/skimage/_shared/interpolation.pyx @@ -7,7 +7,7 @@ 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): + double cval): """Nearest neighbour interpolation at a given position in the image. Parameters @@ -31,7 +31,7 @@ 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): """Bilinear interpolation at a given position in the image. Parameters @@ -65,7 +65,7 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols, cdef inline double get_pixel(double* image, int rows, int cols, int r, int c, - char mode, double cval=0): + char mode, double cval): """Get a pixel from the image, taking wrapping mode into consideration. Parameters diff --git a/skimage/feature/_texture.pyx b/skimage/feature/_texture.pyx index 5fb53e90..70a446bb 100644 --- a/skimage/feature/_texture.pyx +++ b/skimage/feature/_texture.pyx @@ -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(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: diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx index 2851a728..5903380e 100644 --- a/skimage/transform/_warps_cy.pyx +++ b/skimage/transform/_warps_cy.pyx @@ -111,14 +111,17 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1, cdef int rows = img.shape[0] cdef int cols = img.shape[1] + cdef double (*interp_func)(double*, int, int, double, double, + char, double) + if order == 0: + interp_func = nearest_neighbour + elif order == 1: + interp_func = bilinear_interpolation + for tfr in range(out_r): for tfc in range(out_c): _matrix_transform(tfc, tfr, M.data, &c, &r) - if order == 0: - out[tfr, tfc] = nearest_neighbour(img.data, rows, - cols, r, c, mode_c, cval) - elif order == 1: - out[tfr, tfc] = bilinear_interpolation(img.data, rows, - cols, r, c, mode_c, cval) + out[tfr, tfc] = interp_func(img.data, rows, cols, r, c, + mode_c, cval) return out From 7a0e0b8f338e8b758a7ed869fd9be86e011b6a0e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 30 Aug 2012 07:37:07 +0200 Subject: [PATCH 15/24] Simplify mode determination --- skimage/transform/_warps_cy.pyx | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx index 5903380e..18b1a5b0 100644 --- a/skimage/transform/_warps_cy.pyx +++ b/skimage/transform/_warps_cy.pyx @@ -87,13 +87,7 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1, if mode not in ('constant', 'wrap', 'reflect'): raise ValueError("Invalid mode specified. Please use " "`constant`, `wrap` or `reflect`.") - cdef char mode_c - if mode == 'constant': - mode_c = ord('C') - elif mode == 'wrap': - mode_c = ord('W') - elif mode == 'reflect': - mode_c = ord('R') + cdef char mode_c = ord(mode[0].upper()) cdef int out_r, out_c if output_shape is None: From abe5dc3cecd901fa896df7b9b74b8c885e665131 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 30 Aug 2012 09:28:03 +0200 Subject: [PATCH 16/24] Add (bi-)cubic interpolation --- skimage/_shared/interpolation.pxd | 13 +++- skimage/_shared/interpolation.pyx | 101 ++++++++++++++++++++++++-- skimage/transform/_geometric.py | 2 +- skimage/transform/_warps_cy.pyx | 9 ++- skimage/transform/tests/test_warps.py | 2 +- 5 files changed, 113 insertions(+), 14 deletions(-) diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd index 3e8e74e9..40df6524 100644 --- a/skimage/_shared/interpolation.pxd +++ b/skimage/_shared/interpolation.pxd @@ -1,12 +1,19 @@ -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) cdef inline double bilinear_interpolation(double* image, int rows, int cols, double r, double c, char mode, double cval) +cdef inline double cubic_interpolation(double x, double[4] f) + +cdef inline double bicubic_interpolation(double* image, int rows, int cols, + double r, double c, char mode, + double cval) + cdef inline double get_pixel(double* image, int rows, int cols, int r, int c, char mode, double cval) diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx index 63f5b70c..bcebd78d 100644 --- a/skimage/_shared/interpolation.pyx +++ b/skimage/_shared/interpolation.pyx @@ -5,16 +5,17 @@ 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): +cdef inline double nearest_neighbour_interpolation(double* image, int rows, + int cols, double r, + double c, char mode, + double cval): """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. @@ -23,6 +24,11 @@ cdef inline double nearest_neighbour(double* image, int rows, int cols, cval : double Constant value to use for constant mode. + Returns + ------- + value : double + Interpolated value. + """ return get_pixel(image, rows, cols, round(r), round(c), @@ -38,7 +44,7 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols, ---------- image : double array Input image. - rows, cols: int + rows, cols : int Shape of image. r, c : int Position at which to interpolate. @@ -47,6 +53,11 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols, 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,6 +75,79 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols, return (1 - dr) * top + dr * bottom +cdef inline double cubic_interpolation(double x, double[4] f): + """Ccubic 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): + """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'} + Wrapping mode. Constant, Wrap or Reflect. + cval : double + Constant value to use for constant mode. + + Returns + ------- + value : double + Interpolated value. + + """ + + cdef int r0 = r + cdef int c0 = c + if r < 0: + r0 -= 1 + if c < 0: + c0 -= 1 + # scale position to range [0, 1] + cdef double xr = (r - r0 + 1) / 3 + cdef double xc = (c - c0 + 1) / 3 + + cdef double fc[4], fr[4] + + cdef int pr, pc + + for pr in range(r0 - 1, r0 + 3): + + # do row-wise cubic interpolation + for pc in range(c0 - 1, c0 + 3): + fc[pc + 1 - c0] = get_pixel(image, rows, cols, pr, pc, mode, cval) + fr[pr + 1 - r0] = cubic_interpolation(xc, fc) + + # do 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): """Get a pixel from the image, taking wrapping mode into consideration. @@ -72,7 +156,7 @@ cdef inline double get_pixel(double* image, int rows, int cols, int r, int c, ---------- image : double array Input image. - rows, cols: int + rows, cols : int Shape of image. r, c : int Position at which to get the pixel. @@ -81,6 +165,11 @@ cdef inline double get_pixel(double* image, int rows, int cols, int r, int c, 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): diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index f11f1b44..6f895fa7 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -834,7 +834,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, # use fast Cython version for specific parameters fast_modes = ('constant', 'reflect', 'wrap') - if order in (0, 1) and mode in fast_modes and not map_args: + if order in (0, 1, 3) and mode in fast_modes and not map_args: matrix = None if isinstance(inverse_map, HOMOGRAPHY_TRANSFORMS): matrix = inverse_map._matrix diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx index 18b1a5b0..37a577f7 100644 --- a/skimage/transform/_warps_cy.pyx +++ b/skimage/transform/_warps_cy.pyx @@ -5,8 +5,9 @@ cimport numpy as np import numpy as np -from skimage._shared.interpolation cimport (nearest_neighbour, - bilinear_interpolation) +from skimage._shared.interpolation cimport (nearest_neighbour_interpolation, + bilinear_interpolation, + bicubic_interpolation) cdef inline void _matrix_transform(double x, double y, double* H, double *x_, @@ -108,9 +109,11 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1, cdef double (*interp_func)(double*, int, int, double, double, char, double) if order == 0: - interp_func = nearest_neighbour + interp_func = nearest_neighbour_interpolation elif order == 1: interp_func = bilinear_interpolation + elif order == 3: + interp_func = bicubic_interpolation for tfr in range(out_r): for tfc in range(out_c): diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py index d5ffff59..41d1a46a 100644 --- a/skimage/transform/tests/test_warps.py +++ b/skimage/transform/tests/test_warps.py @@ -57,7 +57,7 @@ def test_fast_homography(): tform = ProjectiveTransform(H) coords = warp_coords(tform.inverse, (img.shape[0], img.shape[1])) - for order in range(2): + for order in range(4): for mode in ('constant', 'reflect', 'wrap'): p0 = map_coordinates(img, coords, mode=mode, order=order) p1 = warp(img, tform, mode=mode, order=order) From 1592e47e661c0bea7a965af09a7e3e1b3b020d3f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 30 Aug 2012 09:47:00 +0200 Subject: [PATCH 17/24] Apply clipping also to fast cython implementation --- skimage/transform/_geometric.py | 31 +++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 6f895fa7..130cfc2d 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -834,7 +834,9 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, # use fast Cython version for specific parameters fast_modes = ('constant', 'reflect', 'wrap') - if order in (0, 1, 3) and mode in fast_modes and not map_args: + fast_orders = (0, 1, 3) + + if order in fast_orders and mode in fast_modes and not map_args: matrix = None if isinstance(inverse_map, HOMOGRAPHY_TRANSFORMS): matrix = inverse_map._matrix @@ -851,29 +853,30 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, out = np.dstack(dims) if orig_ndim == 2: out = out[..., 0] - return out - if output_shape is None: - output_shape = ishape + else: # use ndimage.map_coordinates - rows, cols = output_shape[:2] + if output_shape is None: + output_shape = ishape - def coord_map(*args): - return inverse_map(*args, **map_args) + rows, cols = output_shape[:2] - coords = warp_coords(coord_map, (rows, cols, bands)) + def coord_map(*args): + return inverse_map(*args, **map_args) - # Prefilter not necessary for order 1 interpolation - prefilter = order > 1 - mapped = ndimage.map_coordinates(image, coords, prefilter=prefilter, - mode=mode, order=order, cval=cval) + 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() From cb870fd069c16685aef0605c9f2a53fda3d948e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 30 Aug 2012 09:53:41 +0200 Subject: [PATCH 18/24] Fix example of estimate_transform --- skimage/transform/_geometric.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 130cfc2d..95af55d3 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -643,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 From 4dfdc7f74f317fcf57c2ea460ba2e0715bc35d44 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 30 Aug 2012 11:57:05 +0200 Subject: [PATCH 19/24] Update doc string of _warp_fast for bicubic interpolation --- skimage/transform/_warps_cy.pyx | 1 + 1 file changed, 1 insertion(+) diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx index 37a577f7..99eb2500 100644 --- a/skimage/transform/_warps_cy.pyx +++ b/skimage/transform/_warps_cy.pyx @@ -72,6 +72,7 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1, Order of interpolation:: * 0: Nearest-neighbour interpolation. * 1: Bilinear interpolation (default). + * 3: Bicubic interpolation. mode : {'constant', 'reflect', 'wrap'} How to handle values outside the image borders. cval : string From 4cd1f8798b8ff0a615894bb1ac8f4dad7e99b5a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 30 Aug 2012 12:08:21 +0200 Subject: [PATCH 20/24] Add nearest mode for positions outside image --- skimage/_shared/interpolation.pyx | 35 +++++++++++++++------------ skimage/transform/_geometric.py | 7 ++---- skimage/transform/_warps_cy.pyx | 4 +-- skimage/transform/tests/test_warps.py | 2 +- 4 files changed, 25 insertions(+), 23 deletions(-) diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx index bcebd78d..fc247530 100644 --- a/skimage/_shared/interpolation.pyx +++ b/skimage/_shared/interpolation.pyx @@ -19,8 +19,8 @@ cdef inline double nearest_neighbour_interpolation(double* image, int rows, Shape of image. r, c : int Position at which to interpolate. - mode : {'C', 'W', 'R'} - Wrapping mode. Constant, Wrap or Reflect. + mode : {'C', 'W', 'R', 'N'} + Wrapping mode. Constant, Wrap, Reflect or Nearest. cval : double Constant value to use for constant mode. @@ -48,8 +48,8 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols, Shape of image. r, c : int Position at which to interpolate. - mode : {'C', 'W', 'R'} - Wrapping mode. Constant, Wrap or Reflect. + mode : {'C', 'W', 'R', 'N'} + Wrapping mode. Constant, Wrap, Reflect or Nearest. cval : double Constant value to use for constant mode. @@ -111,8 +111,8 @@ cdef inline double bicubic_interpolation(double* image, int rows, int cols, Shape of image. r, c : int Position at which to interpolate. - mode : {'C', 'W', 'R'} - Wrapping mode. Constant, Wrap or Reflect. + mode : {'C', 'W', 'R', 'N'} + Wrapping mode. Constant, Wrap, Reflect or Nearest. cval : double Constant value to use for constant mode. @@ -160,8 +160,8 @@ cdef inline double get_pixel(double* image, int rows, int cols, int r, int c, Shape of image. r, c : int Position at which to get the pixel. - mode : {'C', 'W', 'R'} - Wrapping mode. Constant, Wrap or Reflect. + mode : {'C', 'W', 'R', 'N'} + Wrapping mode. Constant, Wrap, Reflect or Nearest. cval : double Constant value to use for constant mode. @@ -190,28 +190,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', 'R'} + mode : {'W', 'R', 'N'} Whether to wrap or reflect the coordinate if it falls outside [0, dim). """ dim = dim - 1 if mode == 'R': # reflect - if (coord < 0): + if coord < 0: # How many times times does the coordinate wrap? - if ((-coord / dim) % 2 != 0): + if (-coord / dim) % 2 != 0: return dim - (-coord % dim) else: return (-coord % dim) - elif (coord > dim): - if ((coord / dim) % 2 != 0): + elif coord > dim: + if (coord / dim) % 2 != 0: return (dim - (coord % dim)) else: return (coord % dim) elif mode == 'W': # wrap - if (coord < 0): + if coord < 0: return (dim - (-coord % dim)) - elif (coord > dim): + elif coord > dim: return (coord % dim) + elif mode == 'N': # nearest + if coord < 0: + return 0 + elif coord > dim: + return dim return coord diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 95af55d3..02b6c161 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -832,11 +832,8 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, ishape = np.array(image.shape) bands = ishape[2] - # use fast Cython version for specific parameters - fast_modes = ('constant', 'reflect', 'wrap') - fast_orders = (0, 1, 3) - - if order in fast_orders and mode in fast_modes and not map_args: + # use fast Cython version for specific interpolation orders + if order in (0, 1, 3) and not map_args: matrix = None if isinstance(inverse_map, HOMOGRAPHY_TRANSFORMS): matrix = inverse_map._matrix diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx index 99eb2500..bedbfa61 100644 --- a/skimage/transform/_warps_cy.pyx +++ b/skimage/transform/_warps_cy.pyx @@ -86,9 +86,9 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1, cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] M = \ np.ascontiguousarray(H) - if mode not in ('constant', 'wrap', 'reflect'): + if mode not in ('constant', 'wrap', 'reflect', 'nearest'): raise ValueError("Invalid mode specified. Please use " - "`constant`, `wrap` or `reflect`.") + "`constant`, `nearest`, `wrap` or `reflect`.") cdef char mode_c = ord(mode[0].upper()) cdef int out_r, out_c diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py index 41d1a46a..65514073 100644 --- a/skimage/transform/tests/test_warps.py +++ b/skimage/transform/tests/test_warps.py @@ -58,7 +58,7 @@ def test_fast_homography(): coords = warp_coords(tform.inverse, (img.shape[0], img.shape[1])) for order in range(4): - for mode in ('constant', 'reflect', 'wrap'): + for mode in ('constant', 'reflect', 'wrap', 'nearest'): p0 = map_coordinates(img, coords, mode=mode, order=order) p1 = warp(img, tform, mode=mode, order=order) From 146d5a3f5b460de4673f1e63fef1a6b21b8b8b1a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 30 Aug 2012 17:13:30 +0200 Subject: [PATCH 21/24] Remove duplicate subtraction --- skimage/_shared/interpolation.pyx | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx index fc247530..ad8aae4b 100644 --- a/skimage/_shared/interpolation.pyx +++ b/skimage/_shared/interpolation.pyx @@ -123,28 +123,27 @@ cdef inline double bicubic_interpolation(double* image, int rows, int cols, """ - cdef int r0 = r - cdef int c0 = c + cdef int r0 = r - 1 + cdef int c0 = c - 1 if r < 0: r0 -= 1 if c < 0: c0 -= 1 # scale position to range [0, 1] - cdef double xr = (r - r0 + 1) / 3 - cdef double xc = (c - c0 + 1) / 3 + cdef double xr = (r - r0) / 3 + cdef double xc = (c - c0) / 3 cdef double fc[4], fr[4] cdef int pr, pc - for pr in range(r0 - 1, r0 + 3): + # 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) - # do row-wise cubic interpolation - for pc in range(c0 - 1, c0 + 3): - fc[pc + 1 - c0] = get_pixel(image, rows, cols, pr, pc, mode, cval) - fr[pr + 1 - r0] = cubic_interpolation(xc, fc) - - # do cubic interpolation for interpolated values of each row + # cubic interpolation for interpolated values of each row return cubic_interpolation(xr, fr) From 15cc7f1779524acd232caf3ad878a268e7f122b7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 30 Aug 2012 18:14:42 +0200 Subject: [PATCH 22/24] Add biquadratic interpolation --- skimage/_shared/interpolation.pxd | 6 ++- skimage/_shared/interpolation.pyx | 74 ++++++++++++++++++++++++++++++- skimage/transform/_geometric.py | 2 +- skimage/transform/_warps_cy.pyx | 4 ++ 4 files changed, 83 insertions(+), 3 deletions(-) diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd index 40df6524..ef880109 100644 --- a/skimage/_shared/interpolation.pxd +++ b/skimage/_shared/interpolation.pxd @@ -8,8 +8,12 @@ cdef inline double bilinear_interpolation(double* image, int rows, int cols, double r, double c, char mode, double cval) -cdef inline double cubic_interpolation(double x, double[4] f) +cdef inline double quadratic_interpolation(double x, double[3] f) +cdef inline double biquadratic_interpolation(double* image, int rows, int cols, + double r, double c, char mode, + double cval) +cdef inline double cubic_interpolation(double x, double[4] f) cdef inline double bicubic_interpolation(double* image, int rows, int cols, double r, double c, char mode, double cval) diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx index ad8aae4b..3150f41b 100644 --- a/skimage/_shared/interpolation.pyx +++ b/skimage/_shared/interpolation.pyx @@ -75,8 +75,80 @@ 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): + """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): + """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 = round(r) + cdef int c0 = 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): - """Ccubic interpolation. + """Cubic interpolation. Parameters ---------- diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index 02b6c161..270198aa 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -833,7 +833,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, bands = ishape[2] # use fast Cython version for specific interpolation orders - if order in (0, 1, 3) and not map_args: + if order in range(4) and not map_args: matrix = None if isinstance(inverse_map, HOMOGRAPHY_TRANSFORMS): matrix = inverse_map._matrix diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx index bedbfa61..ce400ed6 100644 --- a/skimage/transform/_warps_cy.pyx +++ b/skimage/transform/_warps_cy.pyx @@ -7,6 +7,7 @@ cimport numpy as np import numpy as np from skimage._shared.interpolation cimport (nearest_neighbour_interpolation, bilinear_interpolation, + biquadratic_interpolation, bicubic_interpolation) @@ -72,6 +73,7 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1, Order of interpolation:: * 0: Nearest-neighbour interpolation. * 1: Bilinear interpolation (default). + * 2: Biquadratic interpolation (default). * 3: Bicubic interpolation. mode : {'constant', 'reflect', 'wrap'} How to handle values outside the image borders. @@ -113,6 +115,8 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1, 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 From b2036aee5c81a8a7f7caa9d39711c7100c222da8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Thu, 30 Aug 2012 18:41:00 +0200 Subject: [PATCH 23/24] Add image resize function --- skimage/transform/__init__.py | 2 +- skimage/transform/_warps.py | 51 ++++++++++++++++++++++++++- skimage/transform/tests/test_warps.py | 11 +++++- 3 files changed, 61 insertions(+), 3 deletions(-) diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index 511c8738..0907544b 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -5,4 +5,4 @@ from .integral import * from ._geometric import (warp, warp_coords, estimate_transform, SimilarityTransform, AffineTransform, ProjectiveTransform, PolynomialTransform) -from ._warps import rotate, swirl, homography +from ._warps import resize, rotate, swirl, homography diff --git a/skimage/transform/_warps.py b/skimage/transform/_warps.py index 4e5a7245..1ed0bed9 100644 --- a/skimage/transform/_warps.py +++ b/skimage/transform/_warps.py @@ -1,5 +1,54 @@ import numpy as np -from ._geometric import warp, SimilarityTransform +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.): diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py index 65514073..ac9272c6 100644 --- a/skimage/transform/tests/test_warps.py +++ b/skimage/transform/tests/test_warps.py @@ -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, rotate, +from skimage.transform import (warp, warp_coords, rotate, resize, AffineTransform, ProjectiveTransform, SimilarityTransform) @@ -81,6 +81,15 @@ def test_rotate(): 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()) From b2e4fd6f32630aadf51d1004e34abf49f67aea0d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Fri, 31 Aug 2012 23:46:23 +0200 Subject: [PATCH 24/24] Add parallel execution support --- skimage/_shared/interpolation.pxd | 16 ++++++++-------- skimage/_shared/interpolation.pyx | 16 ++++++++-------- skimage/transform/_warps_cy.pyx | 16 ++++++++++------ skimage/transform/setup.py | 4 +++- 4 files changed, 29 insertions(+), 23 deletions(-) diff --git a/skimage/_shared/interpolation.pxd b/skimage/_shared/interpolation.pxd index ef880109..6038a853 100644 --- a/skimage/_shared/interpolation.pxd +++ b/skimage/_shared/interpolation.pxd @@ -2,23 +2,23 @@ cdef inline double nearest_neighbour_interpolation(double* image, int rows, int cols, double r, double c, char mode, - double cval) + 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) +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) + double cval) nogil -cdef inline double cubic_interpolation(double x, double[4] f) +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) + 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 diff --git a/skimage/_shared/interpolation.pyx b/skimage/_shared/interpolation.pyx index 3150f41b..231ca045 100644 --- a/skimage/_shared/interpolation.pyx +++ b/skimage/_shared/interpolation.pyx @@ -8,7 +8,7 @@ from libc.math cimport ceil, floor, round cdef inline double nearest_neighbour_interpolation(double* image, int rows, int cols, double r, double c, char mode, - double cval): + double cval) nogil: """Nearest neighbour interpolation at a given position in the image. Parameters @@ -37,7 +37,7 @@ cdef inline double nearest_neighbour_interpolation(double* image, int rows, cdef inline double bilinear_interpolation(double* image, int rows, int cols, double r, double c, char mode, - double cval): + double cval) nogil: """Bilinear interpolation at a given position in the image. Parameters @@ -75,7 +75,7 @@ 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): +cdef inline double quadratic_interpolation(double x, double[3] f) nogil: """Quadratic interpolation. Parameters @@ -96,7 +96,7 @@ cdef inline double quadratic_interpolation(double x, double[3] f): cdef inline double biquadratic_interpolation(double* image, int rows, int cols, double r, double c, char mode, - double cval): + double cval) nogil: """Biquadratic interpolation at a given position in the image. Parameters @@ -147,7 +147,7 @@ cdef inline double biquadratic_interpolation(double* image, int rows, int cols, 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 @@ -172,7 +172,7 @@ cdef inline double cubic_interpolation(double x, double[4] f): cdef inline double bicubic_interpolation(double* image, int rows, int cols, double r, double c, char mode, - double cval): + double cval) nogil: """Bicubic interpolation at a given position in the image. Parameters @@ -220,7 +220,7 @@ cdef inline double bicubic_interpolation(double* image, int rows, int cols, cdef inline double get_pixel(double* image, int rows, int cols, int r, int c, - char mode, double cval): + char mode, double cval) nogil: """Get a pixel from the image, taking wrapping mode into consideration. Parameters @@ -251,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. diff --git a/skimage/transform/_warps_cy.pyx b/skimage/transform/_warps_cy.pyx index ce400ed6..9d3dca70 100644 --- a/skimage/transform/_warps_cy.pyx +++ b/skimage/transform/_warps_cy.pyx @@ -5,6 +5,7 @@ cimport numpy as np import numpy as np +from cython.parallel import prange from skimage._shared.interpolation cimport (nearest_neighbour_interpolation, bilinear_interpolation, biquadratic_interpolation, @@ -12,7 +13,7 @@ from skimage._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 @@ -101,8 +102,9 @@ def _warp_fast(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 = out.data cdef int tfr, tfc cdef double r, c @@ -110,7 +112,7 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1, cdef int cols = img.shape[1] cdef double (*interp_func)(double*, int, int, double, double, - char, double) + char, double) nogil if order == 0: interp_func = nearest_neighbour_interpolation elif order == 1: @@ -120,10 +122,12 @@ def _warp_fast(np.ndarray image, np.ndarray H, output_shape=None, int order=1, elif order == 3: interp_func = bicubic_interpolation - for tfr in range(out_r): + 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, M.data, &c, &r) - out[tfr, tfc] = interp_func(img.data, rows, cols, r, c, - mode_c, cval) + out_data[tfr * out_r + tfc] = interp_func(img.data, rows, + cols, r, c, mode_c, cval) return out diff --git a/skimage/transform/setup.py b/skimage/transform/setup.py index 0e415bad..cd2a0f8a 100644 --- a/skimage/transform/setup.py +++ b/skimage/transform/setup.py @@ -20,7 +20,9 @@ def configuration(parent_package='', top_path=None): include_dirs=[get_numpy_include_dirs()]) config.add_extension('_warps_cy', sources=['_warps_cy.c'], - include_dirs=[get_numpy_include_dirs(), '../_shared']) + include_dirs=[get_numpy_include_dirs(), '../_shared'], + extra_compile_args=['-fopenmp'], + extra_link_args=['-fopenmp']) return config