Repository: spacetx/starfish

Issue:
Add support for other types of projections on ImageStack
# Proposal
I think it would be useful to expose other methods of projecting along axes for `ImageStack` (e.g., `DataArray.sum()` and `DataArray.mean()`). This would add support the ExFISH case for projecting across (R, C) before spot segmentation (#1301).

# Implementation
Some possibilities:
1. I think the thing that is the least work is to make a filter, `Filter.Project` that takes projection type as an input argument. This could replace `Filter.MaxProject`.
2. We could do the same as (1), but make it a method of `ImageStack`. I think a lot of the notebooks would have to be updated, as they use `ImageStack.max_proj()`.
3. We could also make one filter component per projection type, but this seems too granular.

As a related aside: I think there should be only one way to max project (i.e., `ImageStack.max_proj()` and `Filter.MaxProject` are redundant).

Thoughts, @dganguli , @ambrosejcarr ?

Discussion:
codecov-io: # [Codecov](https://codecov.io/gh/spacetx/starfish/pull/1342?src=pr&el=h1) Report
> Merging [#1342](https://codecov.io/gh/spacetx/starfish/pull/1342?src=pr&el=desc) into [master](https://codecov.io/gh/spacetx/starfish/commit/f499bc39584c14210aa9489b0a3b0e2650d58de2?src=pr&el=desc) will **decrease** coverage by `20.49%`.
> The diff coverage is `45.94%`.

[![Impacted file tree graph](https://codecov.io/gh/spacetx/starfish/pull/1342/graphs/tree.svg?width=650&token=kliQHejt6T&height=150&src=pr)](https://codecov.io/gh/spacetx/starfish/pull/1342?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master    #1342      +/-   ##
==========================================
- Coverage   89.04%   68.54%   -20.5%     
==========================================
  Files         147      148       +1     
  Lines        5295     5332      +37     
==========================================
- Hits         4715     3655    -1060     
- Misses        580     1677    +1097
```


| [Impacted Files](https://codecov.io/gh/spacetx/starfish/pull/1342?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [starfish/core/image/\_filter/reduce.py](https://codecov.io/gh/spacetx/starfish/pull/1342/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9pbWFnZS9fZmlsdGVyL3JlZHVjZS5weQ==) | `45.94% <45.94%> (ø)` | |
| [starfish/core/experiment/builder/inplace.py](https://codecov.io/gh/spacetx/starfish/pull/1342/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9leHBlcmltZW50L2J1aWxkZXIvaW5wbGFjZS5weQ==) | `0% <0%> (-87.88%)` | :arrow_down: |
| [starfish/test/full\_pipelines/api/test\_dartfish.py](https://codecov.io/gh/spacetx/starfish/pull/1342/diff?src=pr&el=tree#diff-c3RhcmZpc2gvdGVzdC9mdWxsX3BpcGVsaW5lcy9hcGkvdGVzdF9kYXJ0ZmlzaC5weQ==) | `20% <0%> (-80%)` | :arrow_down: |
| [starfish/test/full\_pipelines/api/test\_iss\_api.py](https://codecov.io/gh/spacetx/starfish/pull/1342/diff?src=pr&el=tree#diff-c3RhcmZpc2gvdGVzdC9mdWxsX3BpcGVsaW5lcy9hcGkvdGVzdF9pc3NfYXBpLnB5) | `23.52% <0%> (-76.48%)` | :arrow_down: |
| [starfish/test/full\_pipelines/api/test\_merfish.py](https://codecov.io/gh/spacetx/starfish/pull/1342/diff?src=pr&el=tree#diff-c3RhcmZpc2gvdGVzdC9mdWxsX3BpcGVsaW5lcy9hcGkvdGVzdF9tZXJmaXNoLnB5) | `24.39% <0%> (-75.61%)` | :arrow_down: |
| [starfish/core/experiment/builder/\_\_init\_\_.py](https://codecov.io/gh/spacetx/starfish/pull/1342/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9leHBlcmltZW50L2J1aWxkZXIvX19pbml0X18ucHk=) | `30.35% <0%> (-69.65%)` | :arrow_down: |
| [starfish/core/intensity\_table/overlap.py](https://codecov.io/gh/spacetx/starfish/pull/1342/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9pbnRlbnNpdHlfdGFibGUvb3ZlcmxhcC5weQ==) | `31.37% <0%> (-64.71%)` | :arrow_down: |
| [starfish/core/spacetx\_format/test\_fuzz.py](https://codecov.io/gh/spacetx/starfish/pull/1342/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9zcGFjZXR4X2Zvcm1hdC90ZXN0X2Z1enoucHk=) | `37.5% <0%> (-62.5%)` | :arrow_down: |
| [.../spots/\_detect\_spots/local\_search\_blob\_detector.py](https://codecov.io/gh/spacetx/starfish/pull/1342/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9zcG90cy9fZGV0ZWN0X3Nwb3RzL2xvY2FsX3NlYXJjaF9ibG9iX2RldGVjdG9yLnB5) | `32.2% <0%> (-61.87%)` | :arrow_down: |
| [starfish/core/spacetx\_format/util.py](https://codecov.io/gh/spacetx/starfish/pull/1342/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9zcGFjZXR4X2Zvcm1hdC91dGlsLnB5) | `37.5% <0%> (-58.34%)` | :arrow_down: |
| ... and [65 more](https://codecov.io/gh/spacetx/starfish/pull/1342/diff?src=pr&el=tree-more) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/spacetx/starfish/pull/1342?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/spacetx/starfish/pull/1342?src=pr&el=footer). Last update [f499bc3...b387e52](https://codecov.io/gh/spacetx/starfish/pull/1342?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

kevinyamauchi: @dganguli , I think @ttung suggested that we keep `Filter.max_proj()` and `ImageStack.max_proj()` for a bit for API stability. I would be in favor of removing `ImageStack.max_proj()` and refactoring `Filter.max_proj()` to use `Filter.reduce()`. I think it's worth keeping `Filter.max_proj()` around because maximum intensity projection is such a common operation. Thoughts?

Patch:
diff --git a/docs/source/_static/data_formatting_examples/format_merfish_U2OS_data.py b/docs/source/_static/data_formatting_examples/format_merfish_U2OS_data.py
index f4866fe6..7ff039da 100644
--- a/docs/source/_static/data_formatting_examples/format_merfish_U2OS_data.py
+++ b/docs/source/_static/data_formatting_examples/format_merfish_U2OS_data.py
@@ -6,16 +6,18 @@ The following script formats MERFISH data acquired from cultured U2-OS cells.
 """
 import argparse
 import functools
+import json
 import os
 from typing import IO, Mapping, Tuple, Union
 
 import numpy as np
+import pandas as pd
 from skimage.io import imread
 from slicedimage import ImageFormat
 
+from starfish.core.util.argparse import FsExistsType
 from starfish.experiment.builder import FetchedTile, TileFetcher, write_experiment_json
 from starfish.types import Axes, Coordinates, Number
-from starfish.util.argparse import FsExistsType
 
 SHAPE = {Axes.Y: 2048, Axes.X: 2048}
 
@@ -29,7 +31,7 @@ def cached_read_fn(file_path):
 
 
 class MERFISHTile(FetchedTile):
-    def __init__(self, file_path, r, ch):
+    def __init__(self, file_path, r, ch, coordinates):
         self.file_path = file_path
         # how to index tiles from indices into multi-page tiff
         # key is a tuple of round, chan. val is the index
@@ -51,6 +53,7 @@ class MERFISHTile(FetchedTile):
                     (7, 1): 14}
         self.r = r
         self.ch = ch
+        self._coordinates = coordinates
 
     @property
     def shape(self) -> Mapping[Axes, int]:
@@ -58,21 +61,17 @@ class MERFISHTile(FetchedTile):
 
     @property
     def coordinates(self) -> Mapping[Union[str, Coordinates], Union[Number, Tuple[Number, Number]]]:
-        # FIXME: (dganguli) please provide proper coordinates here.
-        return {
-            Coordinates.X: (0.0, 0.0001),
-            Coordinates.Y: (0.0, 0.0001),
-            Coordinates.Z: (0.0, 0.0001),
-        }
+        return self._coordinates
 
     def tile_data(self) -> IO:
         return cached_read_fn(self.file_path)[self.map[(self.r, self.ch)], :, :]
 
 
 class MERFISHAuxTile(FetchedTile):
-    def __init__(self, file_path):
+    def __init__(self, file_path, coordinates):
         self.file_path = file_path
         self.dapi_index = 17
+        self._coordinates = coordinates
 
     @property
     def shape(self) -> Mapping[Axes, int]:
@@ -80,12 +79,7 @@ class MERFISHAuxTile(FetchedTile):
 
     @property
     def coordinates(self) -> Mapping[Union[str, Coordinates], Union[Number, Tuple[Number, Number]]]:
-        # FIXME: (dganguli) please provide proper coordinates here.
-        return {
-            Coordinates.X: (0.0, 0.0001),
-            Coordinates.Y: (0.0, 0.0001),
-            Coordinates.Z: (0.0, 0.0001),
-        }
+        return self._coordinates
 
     def tile_data(self) -> np.ndarray:
         return cached_read_fn(self.file_path)[self.dapi_index, :, :]
@@ -95,19 +89,48 @@ class MERFISHTileFetcher(TileFetcher):
     def __init__(self, input_dir, is_dapi):
         self.input_dir = input_dir
         self.is_dapi = is_dapi
+        self.coordinates = self.parse_coordinates(input_dir)
+
+    @staticmethod
+    def parse_coordinates(input_dir):
+        filename = os.path.join(input_dir, "stagePos.csv")
+        data = pd.read_csv(filename, names=['y_min', 'x_min'])
+        data['x_max'] = data['x_min'] + 200
+        data['y_max'] = data['y_min'] + 200
+        return data
+
+    def make_coordinates(self, fov):
+        return {
+            Coordinates.X: (
+                float(self.coordinates.loc[fov, 'x_min']),
+                float(self.coordinates.loc[fov, 'x_max'])
+            ),
+            Coordinates.Y: (
+                float(self.coordinates.loc[fov, 'y_min']),
+                float(self.coordinates.loc[fov, 'y_max'])
+            ),
+            Coordinates.Z: (0.0, 0.001)
+        }
 
     def get_tile(self, fov: int, r: int, ch: int, z: int) -> FetchedTile:
         filename = os.path.join(self.input_dir, 'fov_{}.tif'.format(fov))
         file_path = os.path.join(self.input_dir, filename)
         if self.is_dapi:
-            return MERFISHAuxTile(file_path)
+            return MERFISHAuxTile(file_path, self.make_coordinates(fov))
         else:
-            return MERFISHTile(file_path, r, ch)
+            return MERFISHTile(file_path, r, ch, self.make_coordinates(fov))
 
 
 def format_data(input_dir, output_dir):
-    def add_codebook(experiment_json_doc):
-        experiment_json_doc['codebook'] = "codebook.json"
+
+    input_dir = os.path.abspath(input_dir)
+    output_dir = os.path.abspath(output_dir)
+
+    def add_scale_factors(experiment_json_doc):
+        filename = os.path.join(input_dir, "scale_factors.json")
+        with open(filename, 'r') as f:
+            data = json.load(f)
+        experiment_json_doc['extras'] = {"scale_factors": data}
         return experiment_json_doc
 
     num_fovs = 496
@@ -135,7 +158,7 @@ def format_data(input_dir, output_dir):
                           aux_tile_fetcher={
                               'nuclei': MERFISHTileFetcher(input_dir, is_dapi=True),
                           },
-                          postprocess_func=add_codebook,
+                          postprocess_func=add_scale_factors,
                           default_shape=SHAPE
                           )
 
diff --git a/starfish/core/codebook/codebook.py b/starfish/core/codebook/codebook.py
index 875730a3..5270d98c 100644
--- a/starfish/core/codebook/codebook.py
+++ b/starfish/core/codebook/codebook.py
@@ -450,7 +450,7 @@ class Codebook(xr.DataArray):
         # if a feature is all zero, the information should be spread across the channel
         n = array.sizes[Axes.CH.value] * array.sizes[Axes.ROUND.value]
         partitioned_intensity = np.linalg.norm(np.full(n, fill_value=1 / n), ord=norm_order) / n
-        array = array.fillna(partitioned_intensity)
+        array.values[np.logical_not(np.isfinite(array.values))] = partitioned_intensity
 
         return array, norm
 
diff --git a/starfish/core/experiment/experiment.py b/starfish/core/experiment/experiment.py
index 7c25f9ef..12e4c8d4 100644
--- a/starfish/core/experiment/experiment.py
+++ b/starfish/core/experiment/experiment.py
@@ -1,10 +1,8 @@
-import copy
 import json
-import pprint
 from typing import (
+    Any,
     Callable,
     Collection,
-    Dict,
     Iterator,
     List,
     MutableMapping,
@@ -29,6 +27,8 @@ from starfish.core.imagestack.parser.crop import CropParameters
 from starfish.core.spacetx_format import validate_sptx
 from .version import MAX_SUPPORTED_VERSION, MIN_SUPPORTED_VERSION
 
+_SINGLETON = object()
+
 
 class FieldOfView:
     """
@@ -37,7 +37,10 @@ class FieldOfView:
 
     All images can be accessed using
     :py:func:`~starfish.experiment.experiment.FieldOfView.get_image`
-    with the name of the image type. The primary image is accessed using the name
+    with the name of the image type and any which rounds/chs/zplanes you want to include
+    (if none are provided all are included) If the resulting tileset has unaligned tiles
+    (their x/y coordinates do not all match) we return a list on ImageStacks where each stack
+    represents an aligned group of tiles. The primary image is accessed using the name
     :py:attr:`~starfish.experiment.experiment.FieldOFView.PRIMARY_IMAGES`.
 
     Notes
@@ -46,11 +49,6 @@ class FieldOfView:
     They can obtain their auxiliary image dictionary from a dictionary of auxiliary image to
     :py:class:`~slicedimage.TileSet`.
 
-    When a FieldOfView is initialized we parse each :py:class:`~slicedimage.TileSet`
-    into sub groups according to their physical coordinates. Tiles with the same physical
-    coordinates are grouped together into aligned tilesets. If the FieldOfView is properly
-    registered there should only be one aligned subgroup.
-
     The decoding of :py:class:`~slicedimage.TileSet` to
     :py:class:`~starfish.imagestack.imagestack.ImageStack`
     does not happen until the image is accessed. ImageStacks can only be initialized with
@@ -79,9 +77,6 @@ class FieldOfView:
     ) -> None:
         self._images: MutableMapping[str, TileSet] = dict()
         self._name = name
-        self.aligned_coordinate_groups: Dict[str, List[CropParameters]] = dict()
-        for name, tileset in image_tilesets.items():
-            self.aligned_coordinate_groups[name] = CropParameters.parse_coordinate_groups(tileset)
         self._images = image_tilesets
 
     def __repr__(self):
@@ -106,56 +101,16 @@ class FieldOfView:
         return set(self._images.keys())
 
     def show_aligned_image_groups(self) -> None:
-        """
-        Describe the aligned subgroups for each Tileset in this FOV
-
-        Examples
-        --------
-        >>> fov.show_aligned_image_groups()
-        {'nuclei': ' Group 0:  <starfish.ImageStack r={0}, ch={0}, z={0}, (y, x)=(190,270)>',
-        'primary': ' Group 0:  <starfish.ImageStack r={0, 1, 2, 3, 4, 5}, ch={0, 1, '
-        '2}, z={0}, (y, x)=(190, 270)>'}
-
-        The example describes a FieldOfView with two Tilesets (primary and nuclei), because
-        all images have the same (x, y) coordinates, each Tileset has a single aligned subgroup:
-        Group 0.
-        """
-        all_groups = dict()
-        for name, groups in self.aligned_coordinate_groups.items():
-            y_size = self._images[name].default_tile_shape[0]
-            x_size = self._images[name].default_tile_shape[1]
-            info = '\n'.join(
-                f" Group {k}: "
-                f" <starfish.ImageStack "
-                f"r={v._permitted_rounds if v._permitted_rounds else 1}, "
-                f"ch={v._permitted_chs if v._permitted_chs else 1}, "
-                f"z={v._permitted_zplanes if v._permitted_zplanes else 1}, "
-                f"(y, x)={y_size, x_size}>"
-                for k, v in enumerate(groups)
-            )
-            all_groups[name] = f'{info}'
-        pprint.pprint(all_groups)
+        raise DeprecationWarning("This method has been deprecated. Aligned groups are now parsed "
+                                 "as a part of, FieldOfView.get_images() and are determined by the "
+                                 "selected axes provided to the method.")
 
     def iterate_image_type(self, image_type: str) -> Iterator[ImageStack]:
-        """
-        Iterate through the aligned subgroups of the given image type.
-        (ex. primary)
-
-        Parameters
-        ----------
-        image_type : str
-            The name of the image type to iterate through.
+        raise DeprecationWarning("This method has been deprecated. Instead use "
+                                 "FieldOfView.get_images(image_type)")
 
-        Returns
-        --------
-        Iterator
-            An iterator for the aligned subgroups of the image type
-
-        """
-        for aligned_group, _ in enumerate(self.aligned_coordinate_groups[image_type]):
-            yield self.get_image(item=image_type, aligned_group=aligned_group)
-
-    def get_image(self, item: str, aligned_group: int = 0,
+    def get_image(self, item: str,
+                  aligned_group: Any = _SINGLETON,
                   rounds: Optional[Collection[int]] = None,
                   chs: Optional[Collection[int]] = None,
                   zplanes: Optional[Collection[int]] = None,
@@ -163,15 +118,13 @@ class FieldOfView:
                   y: Optional[Union[int, slice]] = None,
                   ) -> ImageStack:
         """
-        Load into memory the Imagestack representation of an aligned image group. If crop parameters
-        provided, first crop the TileSet.
+        Load into memory the first Imagestack representation of an aligned image group. If crop
+        parameters provided, first crop the TileSet.
 
         Parameters
         ----------
         item: str
             The name of the tileset ex. 'primary' or 'nuclei'
-        aligned_group: int
-            The aligned subgroup, default 0
         rounds : Optional[Collection[int]]
             The rounds in the original dataset to load into the ImageStack.  If this is not set,
             then all rounds are loaded into the ImageStack.
@@ -194,16 +147,76 @@ class FieldOfView:
             The instantiated image stack
 
         """
+        if aligned_group != _SINGLETON:
+            raise RuntimeError("The parameter 'aligned_group` is no longer accepted. Please use"
+                               "FieldOfView.get_images() and  provide sets of selected axes "
+                               "instead")
+        stack_iterator = self.get_images(item=item, rounds=rounds,
+                                         chs=chs, zplanes=zplanes, x=x, y=y)
+        return next(stack_iterator)
+
+    def get_images(self, item: str,
+                   rounds: Optional[Collection[int]] = None,
+                   chs: Optional[Collection[int]] = None,
+                   zplanes: Optional[Collection[int]] = None,
+                   x: Optional[Union[int, slice]] = None,
+                   y: Optional[Union[int, slice]] = None,
+                   ) -> Iterator[ImageStack]:
+        """
+        Load into memory the an iterator of aligned Imagestacks for the given tileset and selected
+        axes.
+
+        Parameters
+        ----------
+        item: str
+            The name of the tileset ex. 'primary' or 'nuclei'
+        rounds : Optional[Collection[int]]
+            The rounds in the original dataset to load into the ImageStack/s.  If this is not set,
+            then all rounds are loaded into the ImageStack.
+        chs : Optional[Collection[int]]
+            The channels in the original dataset to load into the ImageStac/s.  If this is not set,
+            then all channels are loaded into the ImageStack.
+        zplanes : Optional[Collection[int]]
+            The z-layers in the original dataset to load into the ImageStack/s.  If this is not set,
+            then all z-layers are loaded into the ImageStack.
+        x : Optional[Union[int, slice]]
+            The x-range in the x-y tile that is loaded into the ImageStack/s.  If this is not set,
+            then the entire x-y tile is loaded into the ImageStack.
+        y : Optional[Union[int, slice]]
+            The y-range in the x-y tile that is loaded into the ImageStack/s.  If this is not set,
+            then the entire x-y tile is loaded into the ImageStack.
 
-        # Get the set of permitted r/ch/z in the aligned group
-        aligned_group_crop_params = copy.copy((self.aligned_coordinate_groups[item][aligned_group]))
-        if rounds or chs or zplanes or x or y:
-            crop_params = CropParameters(permitted_rounds=rounds, permitted_chs=chs,
-                                         permitted_zplanes=zplanes, x_slice=x, y_slice=y)
-            # Combine with any extra crop parameters given
-            aligned_group_crop_params.further_crop(crop_params)
-        return ImageStack.from_tileset(self._images[item],
-                                       crop_parameters=aligned_group_crop_params)
+        Returns
+        -------
+        AlignedImageStackIterator
+            The instantiated ImageStack or list of Imagestacks if the parameters given include
+            multiple aligned groups.
+
+        """
+        # Parse the tileset into aligned groups, only include tiles from selected axes.
+        aligned_groups = CropParameters.parse_aligned_groups(self._images[item],
+                                                             rounds=rounds, chs=chs,
+                                                             zplanes=zplanes,
+                                                             x=x, y=y)
+        aligned_stack_iterator = AlignedImageStackIterator(tileset=self._images[item],
+                                                           aligned_groups=aligned_groups)
+        return aligned_stack_iterator
+
+
+class AlignedImageStackIterator(Iterator[ImageStack]):
+    """Iterator class of AlignedImageStacks."""
+    def __init__(self, tileset: TileSet, aligned_groups: List[CropParameters]):
+        self.size = len(aligned_groups)
+        self.aligned_groups = iter(aligned_groups)
+        self.tileset = tileset
+
+    def __len__(self) -> int:
+        return self.size
+
+    def __next__(self) -> ImageStack:
+        aligned_group = next(self.aligned_groups)
+        stack = ImageStack.from_tileset(self.tileset, aligned_group)
+        return stack
 
 
 class Experiment:
diff --git a/starfish/core/image/_filter/reduce.py b/starfish/core/image/_filter/reduce.py
new file mode 100644
index 00000000..21d65f3d
--- /dev/null
+++ b/starfish/core/image/_filter/reduce.py
@@ -0,0 +1,142 @@
+from typing import (
+    Callable,
+    Iterable,
+    MutableMapping,
+    Sequence,
+    Union
+)
+
+import numpy as np
+
+from starfish.core.imagestack.imagestack import ImageStack
+from starfish.core.types import Axes, Clip, Coordinates, Number
+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
+        If a string is provided, it should correspond to a numpy function that
+        matches the form specified above
+        (i.e., function is resolved: func = getattr(np, func)).
+        Some common examples below:
+        amax: maximum intensity projection (applies numpy.amax)
+        max: maximum intensity projection (this is an alias for amax and 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 : 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
+
+    See Also
+    --------
+    starfish.types.Axes
+
+    """
+
+    def __init__(
+        self, dims: Iterable[Union[Axes, str]], func: Union[str, Callable] = 'max',
+        clip_method: 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,
+            *args,
+    ) -> ImageStack:
+        """Performs the dimension reduction with the specifed function
+
+        Parameters
+        ----------
+        stack : ImageStack
+            Stack to be filtered.
+
+        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)
+
+        # Update the physical coordinates
+        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(stack._data.coords[coord.value])]
+            else:
+                physical_coords[coord] = reduced.coords[coord.value]
+        reduced_stack = ImageStack.from_numpy(reduced.values, coordinates=physical_coords)
+
+        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))
diff --git a/starfish/core/imagestack/imagestack.py b/starfish/core/imagestack/imagestack.py
index 269d99d1..f7378f67 100644
--- a/starfish/core/imagestack/imagestack.py
+++ b/starfish/core/imagestack/imagestack.py
@@ -260,7 +260,7 @@ class ImageStack:
         """
         config = StarfishConfig()
         tileset = Reader.parse_doc(url, baseurl, backend_config=config.slicedimage)
-        coordinate_groups = CropParameters.parse_coordinate_groups(tileset)
+        coordinate_groups = CropParameters.parse_aligned_groups(tileset)
         crop_params = coordinate_groups[aligned_group]
         return cls.from_tileset(tileset, crop_parameters=crop_params)
 
diff --git a/starfish/core/imagestack/parser/crop.py b/starfish/core/imagestack/parser/crop.py
index cf7f1af5..06583c12 100644
--- a/starfish/core/imagestack/parser/crop.py
+++ b/starfish/core/imagestack/parser/crop.py
@@ -2,7 +2,7 @@ from collections import OrderedDict
 from typing import Collection, List, Mapping, MutableSequence, Optional, Sequence, Tuple, Union
 
 import numpy as np
-from slicedimage import TileSet
+from slicedimage import Tile, TileSet
 
 from starfish.core.imagestack.parser import TileCollectionData, TileData, TileKey
 from starfish.core.types import Axes, Coordinates, Number
@@ -55,31 +55,6 @@ class CropParameters:
         if axis_type == Axes.ZPLANE and self._permitted_zplanes:
             self._permitted_zplanes.add(permitted_axis)
 
-    def further_crop(self, further_crop_params) -> "CropParameters":
-        if (further_crop_params._permitted_chs
-                and not further_crop_params._permitted_chs.issubset(self._permitted_chs)):
-            raise ValueError("Can only provide further ch crop params within the existing set of"
-                             "ch crop params")
-        else:
-            self._permitted_chs = further_crop_params._permitted_chs
-
-        if (further_crop_params._permitted_rounds
-                and not further_crop_params._permitted_rounds.issubset(self._permitted_rounds)):
-            raise ValueError("Can only provide further round crop params within the existing set of"
-                             "round crop params")
-        else:
-            self._permitted_rounds = further_crop_params._permitted_rounds
-
-        if (further_crop_params._permitted_zplanes
-                and not further_crop_params._permitted_zplanes.issubset(self._permitted_zplanes)):
-            raise ValueError("Can only provide further zplane crop params within the existing set "
-                             "of zplane crop params")
-        else:
-            self._permitted_zplanes = further_crop_params._permitted_zplanes
-        self._x_slice = further_crop_params._x_slice
-        self._y_slice = further_crop_params._y_slice
-        return self
-
     def filter_tilekeys(self, tilekeys: Collection[TileKey]) -> Collection[TileKey]:
         """
         Filters tilekeys for those that should be included in the resulting ImageStack.
@@ -130,36 +105,101 @@ class CropParameters:
         return start, stop
 
     @staticmethod
-    def parse_coordinate_groups(tileset: TileSet) -> List["CropParameters"]:
-        """Takes a tileset and compares the physical coordinates on each tile to
-         create aligned coordinate groups (groups of tiles that have the same physical coordinates)
+    def parse_aligned_groups(tileset: TileSet,
+                             rounds: Optional[Collection[int]] = None,
+                             chs: Optional[Collection[int]] = None,
+                             zplanes: Optional[Collection[int]] = None,
+                             x: Optional[Union[int, slice]] = None,
+                             y: Optional[Union[int, slice]] = None
+                             ) -> List["CropParameters"]:
+
+        """Takes a tileset and any optional selected axes lists compares the physical coordinates on each
+         tile to create aligned coordinate groups (groups of tiles that have the same physical
+         coordinates)
+
+        Parameters
+        ----------
+        tileset: TileSet
+            The TileSet to parse
+        rounds: Optional[Collection[int]]
+            The rounds in the tileset to include in the final aligned groups. If this is not set,
+            then all rounds are included.
+        chs: Optional[Collection[int]]
+            The chs in the tileset to include in the final aligned groups. If this is not set,
+            then all chs are included.
+        zplanes: Optional[Collection[int]]
+            The zplanes in the tileset to include in the final aligned groups. If this is not set,
+            then all zplanes are included.
+        x: Optional[Union[int, slice]]
+            The x-range in the x-y tile to include in the final aligned groups.  If this is not set,
+            then the entire x-y tile is included.
+        y: Optional[Union[int, slice]]
+            The y-range in the x-y tile to include in the final aligned groups.  If this is not set,
+            then the entire x-y tile is included.
 
          Returns
          -------
-         A list of CropParameters. Each entry describes the r/ch/z values of tiles that are aligned
-         (have matching coordinates)
+         List["CropParameters"]
+             A list of CropParameters. Each entry describes the r/ch/z values of tiles that are
+             aligned (have matching coordinates) and are within the selected_axes if provided.
          """
         coord_groups: OrderedDict[tuple, CropParameters] = OrderedDict()
         for tile in tileset.tiles():
-            x_y_coords = (
-                tile.coordinates[Coordinates.X][0], tile.coordinates[Coordinates.X][1],
-                tile.coordinates[Coordinates.Y][0], tile.coordinates[Coordinates.Y][1]
-            )
-            # A tile with this (x, y) has already been seen, add tile's Indices to CropParameters
-            if x_y_coords in coord_groups:
-                crop_params = coord_groups[x_y_coords]
-                crop_params._add_permitted_axes(Axes.CH, tile.indices[Axes.CH])
-                crop_params._add_permitted_axes(Axes.ROUND, tile.indices[Axes.ROUND])
-                if Axes.ZPLANE in tile.indices:
-                    crop_params._add_permitted_axes(Axes.ZPLANE, tile.indices[Axes.ZPLANE])
-            else:
-                coord_groups[x_y_coords] = CropParameters(
-                    permitted_chs=[tile.indices[Axes.CH]],
-                    permitted_rounds=[tile.indices[Axes.ROUND]],
-                    permitted_zplanes=[tile.indices[Axes.ZPLANE]] if Axes.ZPLANE in tile.indices
-                    else None)
+            if CropParameters.tile_in_selected_axes(tile, rounds, chs, zplanes):
+                x_y_coords = (
+                    tile.coordinates[Coordinates.X][0], tile.coordinates[Coordinates.X][1],
+                    tile.coordinates[Coordinates.Y][0], tile.coordinates[Coordinates.Y][1]
+                )
+                # A tile with this (x, y) has already been seen, add tile's indices to
+                # CropParameters
+                if x_y_coords in coord_groups:
+                    crop_params = coord_groups[x_y_coords]
+                    crop_params._add_permitted_axes(Axes.CH, tile.indices[Axes.CH])
+                    crop_params._add_permitted_axes(Axes.ROUND, tile.indices[Axes.ROUND])
+                    if Axes.ZPLANE in tile.indices:
+                        crop_params._add_permitted_axes(Axes.ZPLANE, tile.indices[Axes.ZPLANE])
+                else:
+                    coord_groups[x_y_coords] = CropParameters(
+                        permitted_chs=[tile.indices[Axes.CH]],
+                        permitted_rounds=[tile.indices[Axes.ROUND]],
+                        permitted_zplanes=[tile.indices[Axes.ZPLANE]]
+                        if Axes.ZPLANE in tile.indices else None,
+                        x_slice=x,
+                        y_slice=y)
         return list(coord_groups.values())
 
+    @staticmethod
+    def tile_in_selected_axes(tile: Tile,
+                              rounds: Optional[Collection[int]] = None,
+                              chs: Optional[Collection[int]] = None,
+                              zplanes: Optional[Collection[int]] = None) -> bool:
+        """
+        Return True if a tile belongs in a list of selected axes.
+
+        Parameters
+        ----------
+        tile:
+            The tile in question
+        rounds: Optional[Collection[int]]
+            The allowed rounds.
+        chs: Optional[Collection[int]]
+            The allowed chs.
+        zplanes: Optional[Collection[int]]
+            The allowed zplanes.
+
+        Returns
+        -------
+        Boolean
+            True if tile belongs with selected axes, False if not.
+        """
+        if rounds and tile.indices[Axes.ROUND] not in rounds:
+            return False
+        if chs and tile.indices[Axes.CH] not in chs:
+            return False
+        if zplanes and tile.indices[Axes.ZPLANE] not in zplanes:
+            return False
+        return True
+
     def crop_shape(self, shape: Mapping[Axes, int]) -> Mapping[Axes, int]:
         """
         Given the shape of the original tile, return the shape of the cropped tile.


Test patch:
diff --git a/starfish/core/experiment/test/test_experiment.py b/starfish/core/experiment/test/test_experiment.py
index 865609c1..0febec5b 100644
--- a/starfish/core/experiment/test/test_experiment.py
+++ b/starfish/core/experiment/test/test_experiment.py
@@ -24,8 +24,8 @@ def round_to_z(r: int) -> Tuple[float, float]:
 NUM_ROUND = 5
 NUM_CH = 2
 NUM_Z = 1
-HEIGHT = 10
-WIDTH = 10
+HEIGHT = 100
+WIDTH = 100
 
 
 def get_aligned_tileset():
@@ -111,25 +111,18 @@ def test_crop_experiment():
 def test_fov_aligned_tileset():
     tilesets = {'primary': get_aligned_tileset(), 'nuclei': get_aligned_tileset()}
     fov = FieldOfView("aligned", tilesets)
-    # Assert only one coordinate group for each aligned stack
-    assert len(fov.aligned_coordinate_groups['primary']) == 1
-    assert len(fov.aligned_coordinate_groups['nuclei']) == 1
+    primary_images = fov.get_images('primary')
+    nuclei_images = fov.get_images('nuclei')
+    # Assert that only one ImageStack in list
+    assert len(primary_images) == 1
+    assert len(nuclei_images) == 1
 
 
 def test_fov_un_aligned_tileset():
     tilesets = {'primary': get_un_aligned_tileset(), 'nuclei': get_un_aligned_tileset()}
     fov = FieldOfView("unaligned", tilesets)
+    primary_images = fov.get_images('primary')
+    nuclei_images = fov.get_images('nuclei')
     # Assert that the number of coordinate groups == NUM_ROUNDS
-    assert len(fov.aligned_coordinate_groups['primary']) == NUM_ROUND
-    assert len(fov.aligned_coordinate_groups['nuclei']) == NUM_ROUND
-
-
-def test_fov_partially_aligned_tileset():
-    # Create combo of aligned and unaligned tiles
-    partially_aligned_tileset = get_aligned_tileset()
-    for tile in get_un_aligned_tileset().tiles():
-        partially_aligned_tileset.add_tile(tile)
-    tileset_dict = {'primary': partially_aligned_tileset}
-    fov = FieldOfView("paritally aligned", tileset_dict)
-    # Assert the number of tile groups is num rounds + 1 aligned group
-    assert len(fov.aligned_coordinate_groups['primary']) == NUM_ROUND + 1
+    assert len(primary_images) == NUM_ROUND
+    assert len(nuclei_images) == NUM_ROUND
diff --git a/starfish/core/image/_filter/test/test_api_contract.py b/starfish/core/image/_filter/test/test_api_contract.py
index 3824450e..792349b9 100644
--- a/starfish/core/image/_filter/test/test_api_contract.py
+++ b/starfish/core/image/_filter/test/test_api_contract.py
@@ -25,6 +25,7 @@ import pytest
 from starfish import ImageStack
 from starfish.core.image import Filter
 from starfish.core.image._filter.max_proj import MaxProject
+from starfish.core.image._filter.reduce import Reduce
 
 methods: Mapping[str, Type] = Filter._algorithm_to_class_map()
 
@@ -53,38 +54,42 @@ def test_all_methods_adhere_to_contract(filter_class):
     # assert isinstance(volume_param, bool), \
     #     f'{filter_class} is_volume must be a bool, not {type(volume_param)}'
 
-    # always emits an Image, even if in_place=True and the resulting filter operates in-place
     data = generate_default_data()
-    try:
-        filtered = instance.run(data, in_place=True)
-    except TypeError:
-        raise AssertionError(f'{filter_class} must accept in_place parameter')
-    assert isinstance(filtered, ImageStack)
-    if filter_class is not MaxProject:
-        # Max Proj does not have an in place option, so we need to skip this assertion
+
+    # Max Proj and Reduce don't have an in_place, n_processes, verbose option,
+    # so we need to skip these tests
+    if filter_class not in [MaxProject, Reduce]:
+        # always emits an Image, even if in_place=True and the resulting filter operates in-place
+        try:
+            filtered = instance.run(data, in_place=True)
+        except TypeError:
+            raise AssertionError(f'{filter_class} must accept in_place parameter')
+        assert isinstance(filtered, ImageStack)
         assert data is filtered, \
             f'{filter_class} should return a reference to the input ImageStack when run in_place'
 
-    # operates out of place
-    data = generate_default_data()
-    filtered = instance.run(data, in_place=False)
-    assert data is not filtered, \
-        f'{filter_class} should output a new ImageStack when run out-of-place'
-
-    # accepts n_processes
-    # TODO shanaxel: verify that this causes more than one process to be generated
-    data = generate_default_data()
-    try:
-        instance.run(data, n_processes=1)
-    except TypeError:
-        raise AssertionError(f'{filter_class} must accept n_processes parameter')
-
-    # accepts verbose, and if passed, prints progress
-    data = generate_default_data()
-    try:
-        instance.run(data, verbose=True)
-    except TypeError:
-        raise AssertionError(f'{filter_class} must accept verbose parameter')
+        # operates out of place
+        data = generate_default_data()
+        filtered = instance.run(data, in_place=False)
+        assert data is not filtered, \
+            f'{filter_class} should output a new ImageStack when run out-of-place'
+
+        # accepts n_processes
+        # TODO shanaxel: verify that this causes more than one process to be generated
+        data = generate_default_data()
+        try:
+            instance.run(data, n_processes=1)
+        except TypeError:
+            raise AssertionError(f'{filter_class} must accept n_processes parameter')
+
+        # accepts verbose, and if passed, prints progress
+        data = generate_default_data()
+        try:
+            instance.run(data, verbose=True)
+        except TypeError:
+            raise AssertionError(f'{filter_class} must accept verbose parameter')
+    else:
+        filtered = instance.run(data)
 
     # output is dtype float and within the expected interval of [0, 1]
     assert filtered.xarray.dtype == np.float32, f'{filter_class} must output float32 data'
diff --git a/starfish/core/image/_filter/test/test_reduce.py b/starfish/core/image/_filter/test/test_reduce.py
new file mode 100644
index 00000000..e51775c5
--- /dev/null
+++ b/starfish/core/image/_filter/test/test_reduce.py
@@ -0,0 +1,126 @@
+from collections import OrderedDict
+
+import numpy as np
+import pytest
+import xarray as xr
+
+from starfish import data
+from starfish import ImageStack
+from starfish.core.image._filter.reduce import Reduce
+from starfish.core.imagestack.test.factories import imagestack_with_coords_factory
+from starfish.core.imagestack.test.imagestack_test_utils import verify_physical_coordinates
+from starfish.types import Axes, PhysicalCoordinateTypes
+
+
+X_COORDS = 1, 2
+Y_COORDS = 4, 6
+Z_COORDS = 1, 3
+
+
+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)
+
+    xr.testing.assert_equal(reduced.xarray, expected_result.xarray)
+
+
+def test_max_projection_preserves_coordinates():
+    e = data.ISS(use_test_data=True)
+    nuclei = e.fov().get_image('nuclei')
+
+    red = Reduce(dims=[Axes.ROUND, Axes.CH, Axes.ZPLANE], func='max')
+    nuclei_proj = red.run(nuclei)
+
+    # Since this data already has only 1 round, 1 ch, 1 zplane
+    # let's just assert that the max_proj operation didn't change anything
+    assert nuclei.xarray.equals(nuclei_proj.xarray)
+
+    stack_shape = OrderedDict([(Axes.ROUND, 3), (Axes.CH, 2),
+                               (Axes.ZPLANE, 3), (Axes.Y, 10), (Axes.X, 10)])
+
+    # Create stack with coordinates, verify coords unaffected by max_poj
+    physical_coords = OrderedDict([(PhysicalCoordinateTypes.X_MIN, X_COORDS[0]),
+                                   (PhysicalCoordinateTypes.X_MAX, X_COORDS[1]),
+                                   (PhysicalCoordinateTypes.Y_MIN, Y_COORDS[0]),
+                                   (PhysicalCoordinateTypes.Y_MAX, Y_COORDS[1]),
+                                   (PhysicalCoordinateTypes.Z_MIN, Z_COORDS[0]),
+                                   (PhysicalCoordinateTypes.Z_MAX, Z_COORDS[1])])
+
+    stack = imagestack_with_coords_factory(stack_shape, physical_coords)
+
+    stack_proj = red.run(stack)
+    expected_z = np.average(Z_COORDS)
+    verify_physical_coordinates(stack_proj, X_COORDS, Y_COORDS, expected_z)