diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py index bd93e959..5fe2cbf0 100644 --- a/skimage/transform/_geometric.py +++ b/skimage/transform/_geometric.py @@ -912,7 +912,7 @@ def _stackcopy(a, b): def warp_coords(coord_map, shape, dtype=np.float64): - """Build the source coordinates for the output pixels of an image warp. + """Build the source coordinates for the output of a 2-D image warp. Parameters ---------- @@ -934,8 +934,9 @@ def warp_coords(coord_map, shape, dtype=np.float64): Notes ----- - This is a lower-level routine that produces the source coordinates used by - `warp()`. + + This is a lower-level routine that produces the source coordinates for 2-D + images 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 @@ -946,7 +947,7 @@ def warp_coords(coord_map, shape, dtype=np.float64): Examples -------- - Produce a coordinate map that Shifts an image up and to the right: + Produce a coordinate map that shifts an image up and to the right: >>> from skimage import data >>> from scipy.ndimage import map_coordinates @@ -989,18 +990,45 @@ def warp_coords(coord_map, shape, dtype=np.float64): def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, - mode='constant', cval=0., reverse_map=None): + mode='constant', cval=0., clip=True): """Warp an image according to a given coordinate transformation. Parameters ---------- - image : 2-D or 3-D array + image : ndarray Input image. - inverse_map : transformation object, callable ``xy = f(xy, **kwargs)``, (3, 3) array - Inverse coordinate map. A function that transforms a (N, 2) array of - ``(row, col)`` coordinates in the *output image* into their corresponding - coordinates in the *source image* (e.g. a transformation object or its - inverse). See example section for usage. + inverse_map : transformation object, callable ``cr = f(cr, **kwargs)``, or ndarray + Inverse coordinate map, which transforms coordinates in the output + images into their corresponding coordinates in the input image. + + There are a number of different options to define this map, depending + on the dimensionality of the input image. A 2-D image can have 2 + dimensions for gray-scale images, or 3 dimensions with color + information. + + - For 2-D images, you can directly pass a transformation object, + e.g. `skimage.transform.SimilarityTransform`, or its inverse. + - For 2-D images, you can pass a ``(3, 3)`` homogeneous + transformation matrix, e.g. + `skimage.transform.SimilarityTransform.params`. + - For 2-D images, a function that transforms a ``(M, 2)`` array of + ``(col, row)`` coordinates in the output image to their + corresponding coordinates in the input image. Extra parameters to + the function can be specified through `map_args`. + - For N-D images, you can directly pass an array of coordinates. + The first dimension specifies the coordinates in the input image, + while the subsequent dimensions determine the position in the + output image. E.g. in case of 2-D images, you need to pass an array + of shape ``(2, rows, cols)``, where `rows` and `cols` determine the + shape of the output image, and the first dimension contains the + ``(row, col)`` coordinate in the input image. + See `scipy.ndimage.map_coordinates` for further documentation. + + Note, that a ``(3, 3)`` matrix is interpreted as a homogeneous + transformation matrix, so you cannot interpolate values from a 3-D + input, if the output is of shape ``(3,)``. + + See example section for usage. map_args : dict, optional Keyword arguments passed to `inverse_map`. output_shape : tuple (rows, cols), optional @@ -1009,18 +1037,23 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, and columns need to be specified. order : int, optional The order of interpolation. The order has to be in the range 0-5: - * 0: Nearest-neighbor - * 1: Bi-linear (default) - * 2: Bi-quadratic - * 3: Bi-cubic - * 4: Bi-quartic - * 5: Bi-quintic + - 0: Nearest-neighbor + - 1: Bi-linear (default) + - 2: Bi-quadratic + - 3: Bi-cubic + - 4: Bi-quartic + - 5: Bi-quintic mode : string, optional Points outside the boundaries of the input are filled according to the given mode ('constant', 'nearest', 'reflect' or 'wrap'). cval : float, optional Used in conjunction with mode 'constant', the value outside the image boundaries. + clip : bool, optional + Whether to clip the output to the float range of ``[0, 1]``, or + ``[-1, 1]`` for input images with negative values. This is enabled by + default, since higher order interpolation may produce values outside + the given input range. Notes ----- @@ -1062,89 +1095,129 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1, >>> warped = warp(image, tform.inverse) + For N-D images you can pass a coordinate array, that specifies the + coordinates in the input image for every element in the output image. E.g. + if you want to rescale a 3-D cube, you can do: + + >>> cube_shape = np.array([30, 30, 30]) + >>> cube = np.random.rand(*cube_shape) + + Setup the coordinate array, that defines the scaling: + + >>> scale = 0.1 + >>> output_shape = (scale * cube_shape).astype(int) + >>> coords0, coords1, coords2 = \ + ... np.mgrid[:output_shape[0], :output_shape[1], :output_shape[2]] + >>> coords = np.array([coords0, coords1, coords2]) + + Assume that the cube contains spatial data, where the first array element + center is at coordinate (0.5, 0.5, 0.5) in real space, i.e. we have to + account for this extra offset when scaling the image: + + >>> coords = (coords + 0.5) / scale - 0.5 + >>> warped = warp(cube, coords) + """ - # Backward API compatibility - if reverse_map is not None: - warnings.warn('`reverse_map` parameter is deprecated and replaced by ' - 'the `inverse_map` parameter.') - inverse_map = reverse_map - if image.ndim < 2 or image.ndim > 3: - raise ValueError("Input must have 2 or 3 dimensions.") - - orig_ndim = image.ndim - image = np.atleast_3d(img_as_float(image)) - ishape = np.array(image.shape) - bands = ishape[2] + image = img_as_float(image) + input_shape = np.array(image.shape) if output_shape is None: - output_shape = ishape + output_shape = input_shape else: output_shape = safe_as_int(output_shape) - out = None - # use fast Cython version for specific interpolation orders and input if order in range(4) and not map_args: + # use fast Cython version for specific interpolation orders and input matrix = None - # inverse_map is a transformation matrix as numpy array if isinstance(inverse_map, np.ndarray) and inverse_map.shape == (3, 3): + # inverse_map is a transformation matrix as numpy array matrix = inverse_map - # inverse_map is a homography elif isinstance(inverse_map, HOMOGRAPHY_TRANSFORMS): + # inverse_map is a homography matrix = inverse_map.params - # inverse_map is the inverse of a homography elif (hasattr(inverse_map, '__name__') and inverse_map.__name__ == 'inverse' and get_bound_method_class(inverse_map) \ in HOMOGRAPHY_TRANSFORMS): + # inverse_map is the inverse of a homography matrix = np.linalg.inv(six.get_method_self(inverse_map).params) if matrix is not None: matrix = matrix.astype(np.double) - # 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] + if image.ndim == 2: + out = _warp_fast(image, matrix, + output_shape=output_shape, + order=order, mode=mode, cval=cval) + elif image.ndim == 3: + 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 out is None: # use ndimage.map_coordinates - rows, cols = output_shape[:2] + if out is None: + # use ndimage.map_coordinates - # inverse_map is a transformation matrix as numpy array - if isinstance(inverse_map, np.ndarray) and inverse_map.shape == (3, 3): + if (isinstance(inverse_map, np.ndarray) + and inverse_map.shape == (3, 3)): + # inverse_map is a transformation matrix as numpy array, + # this is only used for order >= 4. inverse_map = ProjectiveTransform(matrix=inverse_map) - def coord_map(*args): - return inverse_map(*args, **map_args) + if isinstance(inverse_map, np.ndarray): + # inverse_map is directly given as coordinates + coords = inverse_map + else: + # inverse_map is given as function, that transforms (N, 2) + # destination coordinates to their corresponding source + # coordinates. This is only supported for 2(+1)-D images. - coords = warp_coords(coord_map, (rows, cols, bands)) + if image.ndim < 2 or image.ndim > 3: + raise ValueError("Only 2-D images (grayscale or color) are " + "supported, when providing a callable " + "`inverse_map`.") + + def coord_map(*args): + return inverse_map(*args, **map_args) + + if len(input_shape) == 3 and len(output_shape) == 2: + # Input image is 2D and has color channel, but output_shape is + # given for 2-D images. Automatically add the color channel + # dimensionality. + output_shape = (output_shape[0], output_shape[1], + input_shape[2]) + + coords = warp_coords(coord_map, output_shape) # Pre-filtering not necessary for order 0, 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(out, 0, 1) + if clip: + # The spline filters sometimes return results outside [0, 1], + # so clip to ensure valid data - if mode == 'constant' and not (0 <= cval <= 1): - clipped[out == cval] = cval + if np.min(image) < 0: + min_val = -1 + else: + min_val = 0 + max_val = 1 - out = clipped + clipped = np.clip(out, min_val, max_val) - if out.ndim == 3 and orig_ndim == 2: - # remove singleton dimension introduced by atleast_3d - return out[..., 0] - else: - return out + if mode == 'constant' and not (0 <= cval <= 1): + clipped[out == cval] = cval + + out = clipped + + return out diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py index 5c210d9b..0197c715 100644 --- a/skimage/transform/tests/test_warps.py +++ b/skimage/transform/tests/test_warps.py @@ -55,6 +55,36 @@ def test_warp_matrix(): outx = warp(x, matrix, order=5) +def test_warp_nd(): + for dim in range(2, 8): + shape = dim * (5,) + + x = np.zeros(shape, dtype=np.double) + x_c = dim * (2,) + x[x_c] = 1 + refx = np.zeros(shape, dtype=np.double) + refx_c = dim * (1,) + refx[refx_c] = 1 + + coord_grid = dim * (slice(0, 5, 1),) + coords = np.array(np.mgrid[coord_grid]) + 1 + + outx = warp(x, coords, order=0, cval=0) + + assert_array_almost_equal(outx, refx) + + +def test_warp_clip(): + x = 2 * np.ones((5, 5), dtype=np.double) + matrix = np.eye(3) + + outx = warp(x, matrix, order=0, clip=False) + assert_array_almost_equal(x, outx) + + outx = warp(x, matrix, order=0, clip=True) + assert_array_almost_equal(x / 2, outx) + + def test_homography(): x = np.zeros((5, 5), dtype=np.double) x[1, 1] = 1