Merge branch 'enh/roi' into 'master'

Enh/roi

See merge request fsl/fslpy!143

CHANGELOG.rst

@@ -6,6 +6,14 @@ order.
 -------------------------
 
 
+Added
+^^^^^
+
+
+* New :mod:`.image.roi` module, for exracting an ROI of an image, or expanding
+  its field-of-view.
+
+
 Changed
 ^^^^^^^
 

doc/fsl.utils.image.roi.rst

+``fsl.utils.image.roi``
+=======================
+
+.. automodule:: fsl.utils.image.roi
+    :members:
+    :undoc-members:
+    :show-inheritance:

doc/fsl.utils.image.rst

@@ -5,6 +5,7 @@
    :hidden:
 
    fsl.utils.image.resample
+   fsl.utils.image.roi
 
 .. automodule:: fsl.utils.image
     :members:

fsl/utils/image/roi.py

+#!/usr/bin/env python
+#
+# roi.py - Extract an ROI of an image.
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+"""This module provides the :func:`roi` function, which can be used to extract
+a region-of-interest from, or expand the field-of-view of, an :class:`.Image`.
+"""
+
+
+import numpy as np
+
+import fsl.data.image      as fslimage
+import fsl.utils.transform as transform
+
+
+def _normaliseBounds(shape, bounds):
+    """Adjust the given ``bounds`` so that it is compatible with the given
+    ``shape``.
+
+    Bounds must be specified for at least three dimensions - if the shape
+    has more than three dimensions, additional bounds are added.
+
+    A :exc:`ValueError` is raised if the provided bounds are invalid.
+
+    :arg bounds: Sequence of ``(lo, hi)`` bounds - see :func:`roi`.
+    :returns:    An adjusted sequence of bounds.
+    """
+
+    bounds = list(bounds)
+
+    if len(bounds) < 3:
+        raise ValueError('')
+
+    if len(bounds) > len(shape):
+        raise ValueError('')
+
+    for b in bounds:
+        if len(b) != 2 or b[0] >= b[1]:
+            raise ValueError('')
+
+    if len(bounds) < len(shape):
+        for s in shape[len(bounds):]:
+            bounds.append((0, s))
+
+    return bounds
+
+
+def roi(image, bounds):
+    """Extract an ROI from the given ``image`` according to the given
+    ``bounds``.
+
+    This function can also be used to pad, or expand the field-of-view, of an
+    image, by passing in negative low bound values, or high bound values which
+    are larger than the image shape. The padded region will contain zeroes.
+
+    :arg image:  :class:`.Image` object
+
+    :arg bounds: Must be a sequence of tuples, containing the low/high bounds
+                 for each voxel dimension, where the low bound is *inclusive*,
+                 and the high bound is *exclusive*. For 4D images, the bounds
+                 for the fourth dimension are optional.
+
+    :returns:    A new :class:`.Image` object containing the region specified
+                 by the ``bounds``.
+    """
+
+    bounds = _normaliseBounds(image.shape, bounds)
+
+    newshape = [hi - lo for lo, hi in bounds]
+    oldslc   = []
+    newslc   = []
+
+    # Figure out how to slice the input image
+    # data array, and the corresponding slice
+    # in the output data array.
+    for (lo, hi), oldlen, newlen in zip(bounds, image.shape, newshape):
+
+        oldlo = max(lo, 0)
+        oldhi = min(hi, oldlen)
+        newlo = max(0, -lo)
+        newhi = newlo + (oldhi - oldlo)
+
+        oldslc.append(slice(oldlo, oldhi))
+        newslc.append(slice(newlo, newhi))
+
+    oldslc = tuple(oldslc)
+    newslc = tuple(newslc)
+
+    # Copy the ROI into the new data array
+    newdata         = np.zeros(newshape, dtype=image.dtype)
+    newdata[newslc] = image.data[oldslc]
+
+    # Create a new affine for the ROI,
+    # with an appropriate offset along
+    # each spatial dimension
+    oldaff = image.voxToWorldMat
+    offset = [lo for lo, hi in bounds[:3]]
+    offset = transform.scaleOffsetXform([1, 1, 1], offset)
+    newaff = transform.concat(oldaff, offset)
+
+    return fslimage.Image(newdata,
+                          xform=newaff,
+                          header=image.header,
+                          name=image.name + '_roi')

tests/test_image_roi.py

+#!/usr/bin/env python
+#
+# test_image_roi.py -
+#
+# Author: Paul McCarthy <pauldmccarthy@gmail.com>
+#
+
+import pytest
+
+import numpy as np
+
+import fsl.data.image as fslimage
+import fsl.utils.image.roi as roi
+
+
+def test_roi():
+
+    # inshape, bounds, expected outshape, expected affine offset
+    tests = [
+        # 3D image, 3D roi
+        ([10, 10, 10],     [(0, 10), (0, 10), (0, 10)], [10, 10, 10],     [0, 0, 0]),
+        ([10, 10, 10],     [(1, 10), (1, 10), (1, 10)], [ 9,  9,  9],     [1, 1, 1]),
+        ([10, 10, 10],     [(1,  9), (1,  9), (1,  9)], [ 8,  8,  8],     [1, 1, 1]),
+        ([10, 10, 10],     [(3,  5), (3,  5), (3,  5)], [ 2,  2,  2],     [3, 3, 3]),
+        ([10, 10, 10],     [(4,  5), (4,  5), (4,  5)], [ 1,  1,  1],     [4, 4, 4]),
+
+        # 4D image, 3D roi
+        ([10, 10, 10, 10], [(0, 10), (0, 10), (0, 10)], [10, 10, 10, 10], [0, 0, 0]),
+        ([10, 10, 10, 10], [(1, 10), (1, 10), (1, 10)], [ 9,  9,  9, 10], [1, 1, 1]),
+        ([10, 10, 10, 10], [(1,  9), (1,  9), (1,  9)], [ 8,  8,  8, 10], [1, 1, 1]),
+        ([10, 10, 10, 10], [(3,  5), (3,  5), (3,  5)], [ 2,  2,  2, 10], [3, 3, 3]),
+        ([10, 10, 10, 10], [(4,  5), (4,  5), (4,  5)], [ 1,  1,  1, 10], [4, 4, 4]),
+
+        # 4D image, 4D roi
+        ([10, 10, 10, 10], [(0, 10), (0, 10), (0, 10), (0, 10)], [10, 10, 10, 10], [0, 0, 0]),
+        ([10, 10, 10, 10], [(1, 10), (1, 10), (1, 10), (1, 10)], [ 9,  9,  9,  9], [1, 1, 1]),
+        ([10, 10, 10, 10], [(1,  9), (1,  9), (1,  9), (1,  9)], [ 8,  8,  8,  8], [1, 1, 1]),
+        ([10, 10, 10, 10], [(3,  5), (3,  5), (3,  5), (3,  5)], [ 2,  2,  2,  2], [3, 3, 3]),
+        ([10, 10, 10, 10], [(4,  5), (4,  5), (4,  5), (4,  5)], [ 1,  1,  1],     [4, 4, 4]),
+
+        # expanding FOV
+        ([10, 10, 10], [(-5, 15), ( 0, 10), ( 0, 10)], [20, 10, 10], [-5,  0,  0]),
+        ([10, 10, 10], [(-5, 15), (-5, 15), ( 0, 10)], [20, 20, 10], [-5, -5,  0]),
+        ([10, 10, 10], [(-5, 15), (-5, 10), (-5, 15)], [20, 15, 20], [-5, -5, -5]),
+        ([10, 10, 10], [(-5, 15), ( 3, 7),  ( 0, 10)], [20,  4, 10], [-5,  3,  0]),
+    ]
+
+    for inshape, bounds, outshape, offset in tests:
+        data  = np.random.randint(1, 10, inshape)
+        image = fslimage.Image(data, xform=np.eye(4))
+
+        result = roi.roi(image, bounds)
+
+        expaff        = np.eye(4)
+        expaff[:3, 3] = offset
+
+        assert np.all(list(result.shape) == list(outshape))
+        assert np.all(np.isclose(result.voxToWorldMat, expaff))
+
+        oldslc = []
+        newslc = []
+
+        for (lo, hi), oldlen in zip(bounds, inshape):
+            oldslc.append(slice(max(lo, 0), min(hi, oldlen)))
+
+        if len(oldslc) < len(inshape):
+            for d in inshape[len(oldslc):]:
+                oldslc.append(slice(0, d))
+
+        for (lo, hi), slc in zip(bounds, oldslc):
+            if lo < 0:    newlo = -lo
+            else:         newlo =  0
+
+            oldlen = slc.stop - slc.start
+
+            newslc.append(slice(newlo, newlo + oldlen))
+
+        if len(newslc) > len(outshape):
+            newslc = newslc[:len(outshape)]
+
+        assert np.all(data[tuple(oldslc)] == result.data[tuple(newslc)])
+
+    # Error on:
+    #  - not enough bounds
+    #  - too many bounds
+    #  - hi >= lo
+    data  = np.random.randint(1, 10, (10, 10, 10))
+    image = fslimage.Image(data, xform=np.eye(4))
+    with pytest.raises(ValueError): roi.roi(image, [(0, 10), (0, 10)])
+    with pytest.raises(ValueError): roi.roi(image, [(0, 10), (0, 10), (0, 10), (0, 10)])
+    with pytest.raises(ValueError): roi.roi(image, [(5, 5),  (0, 10), (0, 10)])
+    with pytest.raises(ValueError): roi.roi(image, [(6, 5),  (0, 10), (0, 10)])
+    with pytest.raises(ValueError): roi.roi(image, [(0, 10), (5,  5), (0, 10)])
+    with pytest.raises(ValueError): roi.roi(image, [(0, 10), (6,  5), (0, 10)])
+    with pytest.raises(ValueError): roi.roi(image, [(0, 10), (0, 10), (5,  5)])
+    with pytest.raises(ValueError): roi.roi(image, [(0, 10), (0, 10), (6,  5)])