Add Filter.Reduce (general dimension reduction for ImageStack)

starfish/core/image/_filter/reduce.py

@@ -0,0 +1,151 @@
+from copy import deepcopy
+from typing import Callable, Iterable, Optional, Union, Sequence
+
+import numpy as np
+
+from starfish.core.imagestack.imagestack import ImageStack
+from starfish.core.types import Axes, Clip, Coordinates
+from starfish.core.util import click
+from starfish.core.util.dtype import preserve_float_range
+from ._base import FilterAlgorithmBase
+
+
+class Reduce(FilterAlgorithmBase):
+    """
+    Reduces the dimensions of the ImageStack by applying a function
+    along one or more axes.
+
+    Parameters
+    ----------
+    dims : Axes
+        one or more Axes to project over
+    func : Union[str, Callable]
+        function to apply across the dimension(s) specified by dims.
+        If a function is provided, it should follow the form specified by
+        DataArray.reduce():
+        http://xarray.pydata.org/en/stable/generated/xarray.DataArray.reduce.html
+
+        The following strings are valid:
+            max: maximum intensity projection (applies numpy.amax)
+            mean: take the mean across the dim(s) (applies numpy.mean)
+            sum: sum across the dim(s) (applies numpy.sum)
+
+    clip_method : Union[str, Clip]
+        (Default Clip.CLIP) Controls the way that data are scaled to retain skimage dtype
+        requirements that float data fall in [0, 1].
+        Clip.CLIP: data above 1 are set to 1, and below 0 are set to 0
+        Clip.SCALE_BY_IMAGE: data above 1 are scaled by the maximum value, with the maximum
+        value calculated over the entire ImageStack
+        Clip.SCALE_BY_CHUNK: data above 1 are scaled by the maximum value, with the maximum
+        value calculated over each slice, where slice shapes are determined by the group_by
+        parameters
+
+    See Also
+    --------
+    starfish.types.Axes
+
+    """
+
+    def __init__(
+        self, dims: Iterable[Union[Axes, str]], func: Union[str, Callable] = 'max',
+        clip_method: Union[str, Clip] = Clip.CLIP
+    ) -> None:
+
+        self.dims = dims
+        self.clip_method = clip_method
+
+        # If the user provided a string, convert to callable
+        if isinstance(func, str):
+            if func == 'max':
+                func = 'amax'
+            func = getattr(np, func)
+        self.func = func
+
+    _DEFAULT_TESTING_PARAMETERS = {"dims": ['r'], "func": 'max'}
+
+    def run(
+            self,
+            stack: ImageStack,
+            in_place: bool = False,
+            verbose: bool = False,
+            n_processes: Optional[int] = None,
+            *args,
+    ) -> ImageStack:
+        """Performs the dimension reduction with the specifed function
+
+        Parameters
+        ----------
+        stack : ImageStack
+            Stack to be filtered.
+        in_place : bool
+            if True, process ImageStack in-place, otherwise return a new stack
+        verbose : bool
+            if True, report on filtering progress (default = False)
+        n_processes : Optional[int]
+            Number of parallel processes to devote to calculating the filter
+
+        Returns
+        -------
+        ImageStack :
+            If in-place is False, return the results of filter as a new stack. Otherwise return the
+            original stack.
+
+        """
+
+        # Apply the reducing function
+        reduced = stack._data.reduce(self.func, dim=[Axes(dim).value for dim in self.dims])
+
+        # Add the reduced dims back and align with the original stack
+        reduced = reduced.expand_dims(tuple(Axes(dim).value for dim in self.dims))
+        reduced = reduced.transpose(*stack.xarray.dims)
+
+        if self.clip_method == Clip.CLIP:
+            reduced = preserve_float_range(reduced, rescale=False)
+        else:
+            reduced = preserve_float_range(reduced, rescale=True)
+
+        physical_coords: MutableMapping[Coordinates, Sequence[Number]] = {}
+        for axis, coord in (
+                (Axes.X, Coordinates.X),
+                (Axes.Y, Coordinates.Y),
+                (Axes.ZPLANE, Coordinates.Z)):
+            if axis in self.dims:
+                # this axis was projected out of existence.
+                assert coord.value not in reduced.coords
+                physical_coords[coord] = [np.average(self._data.coords[coord.value])]
+            else:
+                physical_coords[coord] = reduced.coords[coord.value]
+        reduced_stack = ImageStack.from_numpy(reduced.values, coordinates=physical_coords)
+
+        if in_place:
+            stack._data = reduced_stack._data
+
+            return stack
+        else:
+            return reduced_stack
+
+    @staticmethod
+    @click.command("Reduce")
+    @click.option(
+        "--dims",
+        type=click.Choice(
+            [Axes.ROUND.value, Axes.CH.value, Axes.ZPLANE.value, Axes.X.value, Axes.Y.value]
+        ),
+        multiple=True,
+        help="The dimensions the Imagestack should max project over."
+             "For multiple dimensions add multiple --dims. Ex."
+             "--dims r --dims c")
+    @click.option(
+        "--func",
+        type=click.Choice(["max", "mean", "sum"]),
+        multiple=False,
+        help="The function to apply across dims"
+             "Valid function names: max, mean, sum."
+    )
+    @click.option(
+        "--clip-method", default=Clip.CLIP, type=Clip,
+        help="method to constrain data to [0,1]. options: 'clip', 'scale_by_image', "
+             "'scale_by_chunk'")
+    @click.pass_context
+    def _cli(ctx, dims, func, clip_method):
+        ctx.obj["component"]._cli_run(ctx, Reduce(dims, func, clip_method))

starfish/core/image/_filter/test/test_reduce.py

@@ -0,0 +1,87 @@
+import numpy as np
+import pytest
+import xarray as xr
+
+from starfish import ImageStack
+from starfish.core.image._filter.reduce import Reduce
+from starfish.types import Axes
+
+
+def make_image_stack():
+    '''
+        Make a test ImageStack
+
+    '''
+
+    # Make the test image
+    test = np.ones((2, 4, 1, 2, 2), dtype='float32') * 0.1
+
+    x = [0, 0, 1, 1]
+    y = [0, 1, 0, 1]
+
+    for i in range(4):
+        test[0, i, 0, x[i], y[i]] = 1
+    test[0, 0, 0, 0, 0] = 0.75
+
+    # Make the ImageStack
+    test_stack = ImageStack.from_numpy(test)
+
+    return test_stack
+
+
+def make_expected_image_stack(func):
+    '''
+        Make the expected image stack result
+    '''
+
+    if func == 'max':
+        reduced = np.array(
+            [[[[[0.75, 0.1],
+                [0.1, 0.1]]],
+                [[[0.1, 1],
+                  [0.1, 0.1]]],
+                [[[0.1, 0.1],
+                  [1, 0.1]]],
+                [[[0.1, 0.1],
+                  [0.1, 1]]]]], dtype='float32'
+        )
+    elif func == 'mean':
+        reduced = np.array(
+            [[[[[0.425, 0.1],
+              [0.1, 0.1]]],
+              [[[0.1, 0.55],
+                [0.1, 0.1]]],
+              [[[0.1, 0.1],
+                [0.55, 0.1]]],
+              [[[0.1, 0.1],
+                [0.1, 0.55]]]]], dtype='float32'
+        )
+    elif func == 'sum':
+        reduced = np.array(
+            [[[[[0.85, 0.2],
+              [0.2, 0.2]]],
+              [[[0.2, 1],
+                [0.2, 0.2]]],
+              [[[0.2, 0.2],
+                [1, 0.2]]],
+              [[[0.2, 0.2],
+                [0.2, 1]]]]], dtype='float32'
+        )
+
+    expected_stack = ImageStack.from_numpy(reduced)
+
+    return expected_stack
+
+
+@pytest.mark.parametrize("func", ['max', 'mean', 'sum'])
+def test_image_stack_reduce(func):
+
+    # Get the test stack and expected result
+    test_stack = make_image_stack()
+    expected_result = make_expected_image_stack(func=func)
+
+    # Filter
+    red = Reduce(dims=[Axes.ROUND], func=func)
+    reduced = red.run(test_stack, in_place=False)
+
+    xr.testing.assert_equal(reduced.xarray, expected_result.xarray)