diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py index 4721e7d1..b0eee512 100644 --- a/skimage/transform/__init__.py +++ b/skimage/transform/__init__.py @@ -3,7 +3,7 @@ from .radon_transform import * from .finite_radon_transform import * from ._project import homography as fast_homography from .integral import * -from ._geometric import (warp, estimate_transform, +from ._geometric import (warp, warp_coords, estimate_transform, SimilarityTransform, AffineTransform, ProjectiveTransform, PolynomialTransform) from ._warps import swirl, homography diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index c055247e..cec591a9 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -698,6 +698,83 @@ 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. @@ -721,7 +798,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, mode : string How to handle values outside the image borders. See `scipy.ndimage.map_coordinates` for detail. - cval : string + cval : float Used in conjunction with mode 'constant', the value outside the image boundaries. @@ -753,30 +830,12 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, if output_shape is None: output_shape = ishape - coords = np.empty(np.r_[3, output_shape], dtype=float) - - ## Construct transformed coordinates - rows, cols = output_shape[:2] - # Reshape grid coordinates into a (P, 2) array of (x, y) pairs - tf_coords = np.indices((cols, rows), dtype=float).reshape(2, -1).T + def coord_transform_fn(*args): + return inverse_map(*args, **map_args) - # Map each (x, y) pair to the source image according to - # the user-provided mapping - tf_coords = inverse_map(tf_coords, **map_args) - - # 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, ...]) - - # colour-coordinate mapping - coords[2, ...] = range(bands) + coords = warp_coords(rows, cols, bands, coord_transform_fn) # Prefilter not necessary for order 1 interpolation prefilter = order > 1 @@ -785,4 +844,10 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, # The spline filters sometimes return results outside [0, 1], # so clip to ensure valid data - return np.clip(mapped.squeeze(), 0, 1) + 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/tests/test_warps.py b/skimage/transform/tests/test_warps.py index f35b7f57..ac241722 100644 --- a/skimage/transform/tests/test_warps.py +++ b/skimage/transform/tests/test_warps.py @@ -1,8 +1,11 @@ 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, homography, fast_homography, - SimilarityTransform, ProjectiveTransform) +from skimage.transform import (warp, warp_coords, fast_homography, + AffineTransform, + ProjectiveTransform, + SimilarityTransform) from skimage import transform as tf, data, img_as_float from skimage.color import rgb2gray @@ -25,16 +28,19 @@ def test_homography(): x = np.zeros((5, 5), dtype=np.uint8) x[1, 1] = 255 x = img_as_float(x) - theta = -np.pi/2 - M = np.array([[np.cos(theta),-np.sin(theta),0], - [np.sin(theta), np.cos(theta),4], - [0, 0, 1]]) - x90 = homography(x, M, order=1) + theta = -np.pi / 2 + M = np.array([[np.cos(theta), - np.sin(theta), 0], + [np.sin(theta), np.cos(theta), 4], + [0, 0, 1]]) + + x90 = warp(x, + inverse_map=ProjectiveTransform(M).inverse, + order=1) assert_array_almost_equal(x90, np.rot90(x)) def test_fast_homography(): - img = rgb2gray(data.lena()) + img = rgb2gray(data.lena()).astype(np.uint8) img = img[:, :100] theta = np.deg2rad(30) @@ -74,5 +80,33 @@ def test_swirl(): assert np.mean(np.abs(image - unswirled)) < 0.01 +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) + + +def test_warp_identity(): + lena = img_as_float(rgb2gray(data.lena())) + assert len(lena.shape) == 2 + assert np.allclose(lena, warp(lena, AffineTransform(rotation=0))) + assert not np.allclose(lena, warp(lena, AffineTransform(rotation=0.1))) + rgb_lena = np.transpose(np.asarray([lena, np.zeros_like(lena), lena]), + (1, 2, 0)) + warped_rgb_lena = warp(rgb_lena, AffineTransform(rotation=0.1)) + assert np.allclose(rgb_lena, warp(rgb_lena, AffineTransform(rotation=0))) + assert not np.allclose(rgb_lena, warped_rgb_lena) + # assert no cross-talk between bands + assert np.all(0 == warped_rgb_lena[:, :, 1]) + + +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) + warped_image1 = map_coordinates(image[:, :, 0], coords[:2]) + + if __name__ == "__main__": run_module_suite()