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] 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([