From 15ffd52299360e0d9de89cbc50a8f70178fcbc12 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Wed, 10 Sep 2014 13:42:57 -0400
Subject: [PATCH 01/16] Improve desciption of inverse_map and add option to
 direclty pass coordinates

---
 skimage/transform/_geometric.py | 50 ++++++++++++++++++++++-----------
 1 file changed, 34 insertions(+), 16 deletions(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index bd93e959..fd34d815 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -996,11 +996,26 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
     ----------
     image : 2-D or 3-D array
         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 ``xy = f(xy, **kwargs)``, ndarray
+        Inverse coordinate map, which transforms coordinates in the *output
+        images* into their corresponding coordinates in the *source image*.
+
+        There are a number of different options to define this map:
+
+         - 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 M-D images, a function that transforms a (N, M) coordinates.
+           In case of 2-D images this means a function that transforms a
+           (N, 2) array of ``(x, y)`` coordinates in the *output image* into
+           their corresponding coordinates in the *source image*. Extra
+           parameters to the function can be specified through `map_args`.
+         - For M-D images, you can directly pass an array of coordinates.
+           See `scipy.ndimage.map_coordinates`. Note, that a (3, 3) matrix
+           is interpreted as a homogeneous transformation matrix.
+
+        See example section for usage.
     map_args : dict, optional
         Keyword arguments passed to `inverse_map`.
     output_shape : tuple (rows, cols), optional
@@ -1009,12 +1024,12 @@ 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').
@@ -1120,14 +1135,17 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
     if out is None:  # use ndimage.map_coordinates
         rows, cols = output_shape[:2]
 
-        # inverse_map is a transformation matrix as numpy array
+        # inverse_map is a transformation matrix as numpy array, this is only
+        # used for order >= 4.
         if isinstance(inverse_map, np.ndarray) and inverse_map.shape == (3, 3):
             inverse_map = ProjectiveTransform(matrix=inverse_map)
 
-        def coord_map(*args):
-            return inverse_map(*args, **map_args)
-
-        coords = warp_coords(coord_map, (rows, cols, bands))
+        if isinstance(inverse_map, np.ndarray):
+            coords = inverse_map
+        else:
+            def coord_map(*args):
+                return inverse_map(*args, **map_args)
+            coords = warp_coords(coord_map, (rows, cols, bands))
 
         # Pre-filtering not necessary for order 0, 1 interpolation
         prefilter = order > 1

From ca9a155cf20738d606c0dba99c1f18f83bb5581f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Thu, 11 Sep 2014 10:28:58 -0400
Subject: [PATCH 02/16] Improve inverse_map description

---
 skimage/transform/_geometric.py | 21 +++++++++++++++------
 1 file changed, 15 insertions(+), 6 deletions(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index fd34d815..5e0e7bc0 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -1006,14 +1006,23 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
            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 M-D images, a function that transforms a (N, M) coordinates.
-           In case of 2-D images this means a function that transforms a
-           (N, 2) array of ``(x, y)`` coordinates in the *output image* into
-           their corresponding coordinates in the *source image*. Extra
+         - For M-D images, a function that transforms a (N, M) coordinate
+           matrix in the output image to their corresponding coordinates in
+           the source image, where N is the total number of pixels in the
+           output image. In case of 2-D images this means a function that
+           transforms a (N, 2) array of ``(x, y)`` coordinates. Extra
            parameters to the function can be specified through `map_args`.
          - For M-D images, you can directly pass an array of coordinates.
-           See `scipy.ndimage.map_coordinates`. Note, that a (3, 3) matrix
-           is interpreted as a homogeneous transformation matrix.
+           The first dimension specifies the coordinates in the source image,
+           while the subsequent dimensions determine the position in the
+           output image. 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 source image. 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 `scipy.ndimage.map_coordinates`
+           for further documentation.
 
         See example section for usage.
     map_args : dict, optional

From 1a31968b5453750ec7447dc0d3056111cb0adc36 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Thu, 11 Sep 2014 11:25:22 -0400
Subject: [PATCH 03/16] Make warp function N-D compatible, and add example to
 doc string

---
 skimage/transform/_geometric.py | 111 +++++++++++++++++++-------------
 1 file changed, 65 insertions(+), 46 deletions(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 5e0e7bc0..647adafe 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -994,35 +994,37 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
 
     Parameters
     ----------
-    image : 2-D or 3-D array
+    image : ndarray
         Input image.
     inverse_map : transformation object, callable ``xy = f(xy, **kwargs)``, ndarray
-        Inverse coordinate map, which transforms coordinates in the *output
-        images* into their corresponding coordinates in the *source image*.
+        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:
+        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 M-D images, a function that transforms a (N, M) coordinate
-           matrix in the output image to their corresponding coordinates in
-           the source image, where N is the total number of pixels in the
-           output image. In case of 2-D images this means a function that
-           transforms a (N, 2) array of ``(x, y)`` coordinates. Extra
-           parameters to the function can be specified through `map_args`.
-         - For M-D images, you can directly pass an array of coordinates.
-           The first dimension specifies the coordinates in the source image,
+         - For 2-D images, a function that transforms a ``(M, 2)`` array of
+           ``(x, y)`` 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. In case of 2-D images, you need to pass an array of
-           shape ``(2, rows, cols)``, where `rows` and `cols` determine 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 source image. 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 `scipy.ndimage.map_coordinates`
-           for further documentation.
+           ``(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
@@ -1086,6 +1088,28 @@ 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:
@@ -1093,20 +1117,14 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
                       '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
@@ -1131,30 +1149,35 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
 
         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]
-
         # inverse_map is a transformation matrix as numpy array, this is only
         # used for order >= 4.
-        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 = ProjectiveTransform(matrix=inverse_map)
 
         if isinstance(inverse_map, np.ndarray):
             coords = inverse_map
         else:
+            if image.ndim < 2 or image.ndim > 3:
+                raise ValueError("Input must have 2 or 3 dimensions.")
+
             def coord_map(*args):
                 return inverse_map(*args, **map_args)
-            coords = warp_coords(coord_map, (rows, cols, bands))
+
+            coords = warp_coords(coord_map, output_shape)
 
         # Pre-filtering not necessary for order 0, 1 interpolation
         prefilter = order > 1
@@ -1170,8 +1193,4 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
 
     out = clipped
 
-    if out.ndim == 3 and orig_ndim == 2:
-        # remove singleton dimension introduced by atleast_3d
-        return out[..., 0]
-    else:
-        return out
+    return out

From 9b5b5b8cc5bc06ece3d5e8ffc81238e49f34b128 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Thu, 11 Sep 2014 11:26:07 -0400
Subject: [PATCH 04/16] Remove deprecated parameter of warp

---
 skimage/transform/_geometric.py | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 647adafe..1c397a1e 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -1111,11 +1111,6 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
     >>> 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
 
     image = img_as_float(image)
     input_shape = np.array(image.shape)

From 27d765abd5f54fe700f5d27acb08d8b6f6b3ce2b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Thu, 11 Sep 2014 11:28:06 -0400
Subject: [PATCH 05/16] Highlight fact that warp_coords is only meant for
 2(+1)-D images

---
 skimage/transform/_geometric.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 1c397a1e..cf099c74 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

From 1b49b494ce30f973ca9d420d62debe267f0a7625 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Fri, 12 Sep 2014 16:51:16 -0400
Subject: [PATCH 06/16] Make error message for wrong input dim more explicit

---
 skimage/transform/_geometric.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index cf099c74..8af5588c 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -1168,7 +1168,9 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
             coords = inverse_map
         else:
             if image.ndim < 2 or image.ndim > 3:
-                raise ValueError("Input must have 2 or 3 dimensions.")
+                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)

From f9009410ba1333ff432bd273de5637bab4c6e6d2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Mon, 15 Sep 2014 23:27:03 -0400
Subject: [PATCH 07/16] Add comments to different conditional paths

---
 skimage/transform/_geometric.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 8af5588c..bba323b8 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -1158,15 +1158,20 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
                 out = np.dstack(dims)
 
     if out is None:  # use ndimage.map_coordinates
-        # inverse_map is a transformation matrix as numpy array, this is only
-        # used for order >= 4.
+        # inverse_map is a transformation matrix as numpy array,
+        # this is only used for order >= 4.
         if (isinstance(inverse_map, np.ndarray)
                 and inverse_map.shape == (3, 3)):
             inverse_map = ProjectiveTransform(matrix=inverse_map)
 
         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.
+
             if image.ndim < 2 or image.ndim > 3:
                 raise ValueError("Only 2-D images (grayscale or color) are "
                                  "supported, when providing a callable "
@@ -1179,6 +1184,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
 
         # 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)
 

From 9b4f2d903202f76e70b1d1ff4e790355cc1e3836 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Mon, 15 Sep 2014 23:34:33 -0400
Subject: [PATCH 08/16] Add support for 2(+1)D images and given output_shape
 for 2D images

---
 skimage/transform/_geometric.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index bba323b8..4a1ce6a4 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -1180,6 +1180,13 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
             def coord_map(*args):
                 return inverse_map(*args, **map_args)
 
+            # Input image is 2D and has color channel, but output_shape is
+            # given for 2-D images. Automatically add the color channel
+            # dimensionality.
+            if len(input_shape) == 3 and len(output_shape) == 2:
+                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

From da36750d55e4c8949fec0a64ff705d2c0183ab8c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Tue, 16 Sep 2014 08:19:48 -0400
Subject: [PATCH 09/16] Make row, col ordering clearer in doc string

---
 skimage/transform/_geometric.py | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 4a1ce6a4..20d104e6 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -997,7 +997,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
     ----------
     image : ndarray
         Input image.
-    inverse_map : transformation object, callable ``xy = f(xy, **kwargs)``, ndarray
+    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.
 
@@ -1008,12 +1008,13 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
 
          - 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, 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
-           ``(x, y)`` coordinates in the output image to their corresponding
-           coordinates in the input image. Extra parameters to the function
-           can be specified through `map_args`.
+           ``(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
@@ -1025,7 +1026,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
 
        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, )``.
+       input, if the output is of shape ``(3,)``.
 
         See example section for usage.
     map_args : dict, optional

From 6ff1bf26eea934441a043d0a7065c309e0b99600 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Thu, 18 Sep 2014 20:21:48 -0400
Subject: [PATCH 10/16] Move comments into if scope

---
 skimage/transform/_geometric.py | 22 ++++++++++++----------
 1 file changed, 12 insertions(+), 10 deletions(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 20d104e6..a3dd2a59 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -1124,24 +1124,24 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
 
     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:
@@ -1158,11 +1158,13 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
                                            order=order, mode=mode, cval=cval))
                 out = np.dstack(dims)
 
-    if out is None:  # use ndimage.map_coordinates
-        # inverse_map is a transformation matrix as numpy array,
-        # this is only used for order >= 4.
+    if out is None:
+        # use ndimage.map_coordinates
+
         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)
 
         if isinstance(inverse_map, np.ndarray):
@@ -1181,10 +1183,10 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
             def coord_map(*args):
                 return inverse_map(*args, **map_args)
 
-            # Input image is 2D and has color channel, but output_shape is
-            # given for 2-D images. Automatically add the color channel
-            # dimensionality.
             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])
 

From 6e26f3a7c76d52acac3cfa850a16d9571640d5cd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Thu, 18 Sep 2014 20:27:38 -0400
Subject: [PATCH 11/16] Fix indendation of parameter description

---
 skimage/transform/_geometric.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index a3dd2a59..7813fc35 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -1010,7 +1010,7 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
            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`
+           `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
@@ -1024,9 +1024,9 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
            ``(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,)``.
+        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

From 46305eeeece158366c75bbf7b0ea61045ea72c20 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Thu, 18 Sep 2014 20:36:34 -0400
Subject: [PATCH 12/16] Remove legacy parameter

---
 skimage/transform/_geometric.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 7813fc35..5bc3c5e7 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -990,7 +990,7 @@ 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.):
     """Warp an image according to a given coordinate transformation.
 
     Parameters

From fdf0908ab73a6b8c249e0bc8f765d0f5bd93c2da Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Thu, 18 Sep 2014 20:49:25 -0400
Subject: [PATCH 13/16] Add option to skip clipping of output, and detect if
 input float image has negative values

---
 skimage/transform/_geometric.py       | 26 +++++++++++++++++++-------
 skimage/transform/tests/test_warps.py | 11 +++++++++++
 2 files changed, 30 insertions(+), 7 deletions(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 5bc3c5e7..0264adf0 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -990,7 +990,7 @@ 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.):
+         mode='constant', cval=0., clip=True):
     """Warp an image according to a given coordinate transformation.
 
     Parameters
@@ -1049,6 +1049,10 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
     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]``,
+        since higher order interpolation may produce values outside the
+        given input range.
 
     Notes
     -----
@@ -1198,13 +1202,21 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, 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 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..fcdd0a87 100644
--- a/skimage/transform/tests/test_warps.py
+++ b/skimage/transform/tests/test_warps.py
@@ -55,6 +55,17 @@ def test_warp_matrix():
     outx = warp(x, matrix, order=5)
 
 
+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

From f5c627493ebda4dc47ec3c96d61b4995c7f526ac Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Thu, 18 Sep 2014 20:52:55 -0400
Subject: [PATCH 14/16] Extend description of clipping behavior

---
 skimage/transform/_geometric.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/skimage/transform/_geometric.py b/skimage/transform/_geometric.py
index 0264adf0..5fe2cbf0 100644
--- a/skimage/transform/_geometric.py
+++ b/skimage/transform/_geometric.py
@@ -1050,9 +1050,10 @@ def warp(image, inverse_map=None, map_args={}, output_shape=None, order=1,
         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]``,
-        since higher order interpolation may produce values outside the
-        given input range.
+        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
     -----

From 877a7ba1099e6aa0e2e490cd941d3b61d9d6a83a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Thu, 18 Sep 2014 21:03:56 -0400
Subject: [PATCH 15/16] Add test for N-D warping

---
 skimage/transform/tests/test_warps.py | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py
index fcdd0a87..da4f0beb 100644
--- a/skimage/transform/tests/test_warps.py
+++ b/skimage/transform/tests/test_warps.py
@@ -55,6 +55,25 @@ 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 * (np.arange(5),)
+        coords = np.array(np.meshgrid(*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)

From d367ff479b745c395f1ab130e97ea5ae35cae640 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= <jsch@demuc.de>
Date: Fri, 19 Sep 2014 13:34:00 -0400
Subject: [PATCH 16/16] Use for mgrid for Py26 compatibility

---
 skimage/transform/tests/test_warps.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/skimage/transform/tests/test_warps.py b/skimage/transform/tests/test_warps.py
index da4f0beb..0197c715 100644
--- a/skimage/transform/tests/test_warps.py
+++ b/skimage/transform/tests/test_warps.py
@@ -66,8 +66,8 @@ def test_warp_nd():
         refx_c = dim * (1,)
         refx[refx_c] = 1
 
-        coord_grid = dim * (np.arange(5),)
-        coords = np.array(np.meshgrid(*coord_grid)) + 1
+        coord_grid = dim * (slice(0, 5, 1),)
+        coords = np.array(np.mgrid[coord_grid]) + 1
 
         outx = warp(x, coords, order=0, cval=0)