Merge pull request #511 from ankit-maverick/resample

Resampling of nD arrays

skimage/measure/__init__.py

@@ -3,6 +3,7 @@ from ._regionprops import regionprops, perimeter
 from ._structural_similarity import structural_similarity
 from ._polygon import approximate_polygon, subdivide_polygon
 from .fit import LineModel, CircleModel, EllipseModel, ransac
+from ._sum_blocks import sum_blocks
 
 
 __all__ = ['find_contours',
@@ -14,4 +15,5 @@ __all__ = ['find_contours',
            'LineModel',
            'CircleModel',
            'EllipseModel',
-           'ransac']
+           'ransac',
+           'sum_blocks']

skimage/measure/_sum_blocks.py

@@ -0,0 +1,38 @@
+def sum_blocks(array, factors):
+    """Sums the elements in blocks of integer factors and pads the original
+    array with zeroes if the dimensions are not perfectly divisible by factors.
+
+    This function is different from resize and rescale in transform._warps in
+    the sense that they use interpolation to upsample or downsample on a 2D
+    array, while this function performs only dawnsampling but on any
+    n-dimensional array and returns the sum of elements in a block of size
+    factors in the original array.
+
+    Parameters
+    ----------
+    array : ndarray
+        Input n-dimensional array.
+    factors: tuple
+        Tuple containing integer values representing block length along each
+        axis.
+
+    Returns
+    -------
+    array : ndarray
+        Downsampled array with same number of dimensions as that of input
+        array.
+
+    Example
+    -------
+    >>> a = np.arange(15).reshape(3, 5)
+    >>> a
+    array([[ 0,  1,  2,  3,  4],
+           [ 5,  6,  7,  8,  9],
+           [10, 11, 12, 13, 14]])
+    >>> sum_blocks(a, (2,3))
+    array([[21, 24],
+           [33, 27]])
+
+    """
+    from ..transform._warps import _downsample
+    return _downsample(array, factors)

skimage/measure/tests/test_sum_blocks.py

@@ -0,0 +1,16 @@
+import numpy as np
+from numpy.testing import assert_array_equal
+from skimage.measure._sum_blocks import sum_blocks
+
+def test_downsample_sum_blocks():
+    """Verifying downsampling of an array with expected result in sum mode"""
+    image1 = np.arange(4*6).reshape(4, 6)
+    out1 = sum_blocks(image1, (2, 3))
+    expected1 = np.array([[ 24,  42],
+                          [ 96, 114]])
+    assert_array_equal(expected1, out1)
+    image2 = np.arange(5*8).reshape(5, 8)
+    out2 = sum_blocks(image2, (3, 3))
+    expected2 = np.array([[ 81, 108,  87],
+                          [174, 192, 138]])
+    assert_array_equal(expected2, out2)

skimage/transform/__init__.py

@@ -9,7 +9,7 @@ from ._geometric import (warp, warp_coords, estimate_transform,
                          SimilarityTransform, AffineTransform,
                          ProjectiveTransform, PolynomialTransform,
                          PiecewiseAffineTransform)
-from ._warps import swirl, resize, rotate, rescale
+from ._warps import swirl, resize, rotate, rescale, downscale_local_means
 from .pyramids import (pyramid_reduce, pyramid_expand,
                        pyramid_gaussian, pyramid_laplacian)
 
@@ -40,6 +40,7 @@ __all__ = ['hough_circle',
            'resize',
            'rotate',
            'rescale',
+           'downscale_local_means',
            'pyramid_reduce',
            'pyramid_expand',
            'pyramid_gaussian',

skimage/transform/_warps.py

@@ -1,11 +1,19 @@
 import numpy as np
 from scipy import ndimage
+
 from ._geometric import warp, SimilarityTransform, AffineTransform
+from skimage.util.shape import view_as_blocks, _pad_asymmetric_zeros
 
 
 def resize(image, output_shape, order=1, mode='constant', cval=0.):
     """Resize image to match a certain size.
 
+    Resize performs interpolation to upsample or downsample 2D arrays. For
+    downsampling any n-dimensional array by performing arithmetic sum or
+    arithmetic mean, see measure._sum_blocks.sum_blocks and
+    transform._warps.downscale_local_means respectively.
+
+
     Parameters
     ----------
     image : ndarray
@@ -87,6 +95,11 @@ def resize(image, output_shape, order=1, mode='constant', cval=0.):
 def rescale(image, scale, order=1, mode='constant', cval=0.):
     """Scale image by a certain factor.
 
+    Rescale performs interpolation to upsample or downsample 2D arrays. For
+    downsampling any n-dimensional array by performing arithmetic sum or
+    arithmetic mean, see measure._sum_blocks.sum_blocks and
+    transform._warps.downscale_local_means respectively.
+
     Parameters
     ----------
     image : ndarray
@@ -283,3 +296,90 @@ def swirl(image, center=None, strength=1, radius=100, rotation=0,
     return warp(image, _swirl_mapping, map_args=warp_args,
                 output_shape=output_shape,
                 order=order, mode=mode, cval=cval)
+
+
+def _downsample(array, factors, sum=True):
+    """Performs downsampling with integer factors.
+
+    Parameters
+    ----------
+    array : ndarray
+        Input n-dimensional array.
+    factors: tuple
+        Tuple containing downsampling factor along each axis.
+    sum : bool
+        If True, downsampled element is the sum of its corresponding
+        constituent elements in the input array. Default is True.
+
+    Returns
+    -------
+    array : ndarray
+        Downsampled array with same number of dimensions as that of input
+        array.
+
+    """
+
+    pad_size = []
+    if len(factors) != array.ndim:
+        raise ValueError("'factors' must have the same length "
+                         "as 'array.shape'")
+    else:
+        for i in range(len(factors)):
+            if array.shape[i] % factors[i] != 0:
+                pad_size.append(factors[i] - (array.shape[i] % factors[i]))
+            else:
+                pad_size.append(0)
+
+    for i in range(len(pad_size)):
+        array = _pad_asymmetric_zeros(array, pad_size[i], i)
+
+    out = view_as_blocks(array, factors)
+    block_shape = out.shape
+
+    if sum:
+        for i in range(len(block_shape) // 2):
+            out = out.sum(-1)
+    else:
+        for i in range(len(block_shape) // 2):
+            out = out.mean(-1)
+    return out
+
+
+def downscale_local_means(array, factors):
+    """Downsamples the array in blocks of input integer factors after padding
+    the original array with zeroes if the dimensions are not perfectly
+    divisible by factors and replaces it with mean i.e. average value.
+
+    This function is different from resize and rescale in the sense that they
+    use interpolation to upsample or downsample on a 2D array, while this
+    function performs only dawnsampling but on any n-dimensional array and
+    returns the arithmetic mean of elements in a block of size factors in the
+    original array.
+
+    Parameters
+    ----------
+    array : ndarray
+        Input n-dimensional array.
+    factors: tuple
+        Tuple containing integer values representing block length along each
+        axis.
+
+    Returns
+    -------
+    array : ndarray
+        Downsampled array with same number of dimensions as that of input
+        array.
+
+    Example
+    -------
+    >>> a = np.arange(15).reshape(3, 5)
+    >>> a
+    array([[ 0,  1,  2,  3,  4],
+           [ 5,  6,  7,  8,  9],
+           [10, 11, 12, 13, 14]])
+    >>> downscale_local_means(a, (2,3))
+    array([[3.5, 4.],
+           [5.5, 4.5]])
+
+    """
+    return _downsample(array, factors, False)

skimage/transform/tests/test_warps.py

@@ -1,11 +1,12 @@
-from numpy.testing import assert_array_almost_equal, run_module_suite
+from numpy.testing import assert_array_almost_equal, run_module_suite, assert_array_equal
 import numpy as np
 from scipy.ndimage import map_coordinates
 
 from skimage.transform import (warp, warp_coords, rotate, resize, rescale,
                                AffineTransform,
                                ProjectiveTransform,
-                               SimilarityTransform)
+                               SimilarityTransform,
+                               downscale_local_means)
 from skimage import transform as tf, data, img_as_float
 from skimage.color import rgb2gray
 
@@ -194,5 +195,19 @@ def test_warp_coords_example():
     map_coordinates(image[:, :, 0], coords[:2])
 
 
+def test_downscale_local_means():
+    """Verifying downsampling of an array with expected result in mean mode"""
+    image1 = np.arange(4*6).reshape(4, 6)
+    out1 = downscale_local_means(image1, (2, 3))
+    expected1 = np.array([[  4.,   7.],
+                          [ 16.,  19.]])
+    assert_array_equal(expected1, out1)
+    image2 = np.arange(5*8).reshape(5, 8)
+    out2 = downscale_local_means(image2, (4, 5))
+    expected2 = np.array([[ 14. ,  10.8],
+                          [  8.5,   5.7]])
+    assert_array_equal(expected2, out2)
+
+
 if __name__ == "__main__":
     run_module_suite()

skimage/util/shape.py

@@ -230,3 +230,15 @@ def view_as_windows(arr_in, window_shape):
     arr_out = as_strided(arr_in, shape=new_shape, strides=new_strides)
 
     return arr_out
+
+
+def _pad_asymmetric_zeros(arr, pad_amt, axis=-1):
+    """Pads `arr` with zeros by `pad_amt` along specified `axis`"""
+    if axis == -1:
+        axis = arr.ndim - 1
+
+    zeroshape = tuple([x if i != axis else pad_amt
+                       for (i, x) in enumerate(arr.shape)])
+
+    return np.concatenate((arr, np.zeros(zeroshape, dtype=arr.dtype)),
+                          axis=axis)