Repository: spacetx/starfish

Issue:
Error running PxPx decoding when setting normalize=True for ZeroByChannelMagnitude filter
When running the DARTFISH notebook with: 

z_filt = Filter.ZeroByChannelMagnitude(thresh=.05, normalize=True)

When calling the pixel spot decoder, it returns a value error. 

Running Ubuntu 18.04.2 LTS
Python3.6
starfish-0.1.1+1.g30357c1
slicedimage-3.1.0

Code to reproduce error:

``` python-traceback
import numpy as np
import os
from starfish import Experiment
from starfish import data, FieldOfView
from starfish.types import Features, Axes
from starfish import IntensityTable
from starfish.image import Filter
from starfish.spots import DetectPixels

use_test_data = os.getenv("USE_TEST_DATA") is not None
experiment = data.DARTFISH(use_test_data=use_test_data)
imgs = experiment.fov().get_image(FieldOfView.PRIMARY_IMAGES)
sc_filt = Filter.Clip(p_max=100, expand_dynamic_range=True)
norm_imgs = sc_filt.run(imgs)
z_filt = Filter.ZeroByChannelMagnitude(thresh=.05, normalize=True)
filtered_imgs = z_filt.run(norm_imgs)

magnitude_threshold = 0.5
area_threshold = (5, 30)
distance_threshold = 3

psd = DetectPixels.PixelSpotDecoder(
    codebook=experiment.codebook,
    metric='euclidean',
    distance_threshold=distance_threshold,
    magnitude_threshold=magnitude_threshold,
    min_area=area_threshold[0],
    max_area=area_threshold[1]
)

initial_spot_intensities, results = psd.run(filtered_imgs)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-1-f25b84485c47> in <module>
     29 )
     30 
---> 31 initial_spot_intensities, results = psd.run(filtered_imgs)

~/starfish/.sfenv/lib/python3.6/site-packages/starfish/core/pipeline/algorithmbase.py in helper(*args, **kwargs)
     44         @functools.wraps(func)
     45         def helper(*args, **kwargs):
---> 46             result = func(*args, **kwargs)
     47             # Scenario 1, Filtering, ApplyTransform
     48             if isinstance(result, ImageStack):

~/starfish/.sfenv/lib/python3.6/site-packages/starfish/core/spots/_detect_pixels/pixel_spot_decoder.py in run(self, primary_image, n_processes, *args)
     78             min_intensity=self.magnitude_threshold,
     79             norm_order=self.norm_order,
---> 80             metric=self.metric
     81         )
     82         caf = CombineAdjacentFeatures(

~/starfish/.sfenv/lib/python3.6/site-packages/starfish/core/codebook/codebook.py in decode_metric(self, intensities, max_distance, min_intensity, norm_order, metric)
    563 
    564         metric_outputs, targets = self._approximate_nearest_code(
--> 565             norm_codes, norm_intensities, metric=metric)
    566 
    567         # only targets with low distances and high intensities should be retained

~/starfish/.sfenv/lib/python3.6/site-packages/starfish/core/codebook/codebook.py in _approximate_nearest_code(norm_codes, norm_intensities, metric)
    488         # reshape into traces
    489         nn = NearestNeighbors(n_neighbors=1, algorithm='ball_tree', metric=metric).fit(linear_codes)
--> 490         metric_output, indices = nn.kneighbors(linear_features)
    491         gene_ids = np.ravel(norm_codes.indexes[Features.TARGET].values[indices])
    492 

~/starfish/.sfenv/lib/python3.6/site-packages/sklearn/neighbors/base.py in kneighbors(self, X, n_neighbors, return_distance)
    400         if X is not None:
    401             query_is_train = False
--> 402             X = check_array(X, accept_sparse='csr')
    403         else:
    404             query_is_train = True

~/starfish/.sfenv/lib/python3.6/site-packages/sklearn/utils/validation.py in check_array(array, accept_sparse, accept_large_sparse, dtype, order, copy, force_all_finite, ensure_2d, allow_nd, ensure_min_samples, ensure_min_features, warn_on_dtype, estimator)
    540         if force_all_finite:
    541             _assert_all_finite(array,
--> 542                                allow_nan=force_all_finite == 'allow-nan')
    543 
    544     if ensure_min_samples > 0:

~/starfish/.sfenv/lib/python3.6/site-packages/sklearn/utils/validation.py in _assert_all_finite(X, allow_nan)
     54                 not allow_nan and not np.isfinite(X).all()):
     55             type_err = 'infinity' if allow_nan else 'NaN, infinity'
---> 56             raise ValueError(msg_err.format(type_err, X.dtype))
     57     # for object dtype data, we only check for NaNs (GH-13254)
     58     elif X.dtype == np.dtype('object') and not allow_nan:

ValueError: Input contains NaN, infinity or a value too large for dtype('float32').
```

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

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

```diff
@@             Coverage Diff             @@
##           master    #1371       +/-   ##
===========================================
- Coverage   89.19%   68.61%   -20.59%     
===========================================
  Files         148      148               
  Lines        5368     5368               
===========================================
- Hits         4788     3683     -1105     
- Misses        580     1685     +1105
```


| [Impacted Files](https://codecov.io/gh/spacetx/starfish/pull/1371?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [starfish/core/codebook/codebook.py](https://codecov.io/gh/spacetx/starfish/pull/1371/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9jb2RlYm9vay9jb2RlYm9vay5weQ==) | `79.87% <100%> (-16.99%)` | :arrow_down: |
| [starfish/core/experiment/builder/inplace.py](https://codecov.io/gh/spacetx/starfish/pull/1371/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/1371/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/1371/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/1371/diff?src=pr&el=tree#diff-c3RhcmZpc2gvdGVzdC9mdWxsX3BpcGVsaW5lcy9hcGkvdGVzdF9tZXJmaXNoLnB5) | `24.39% <0%> (-75.61%)` | :arrow_down: |
| [...fish/core/imagestack/test/factories/all\_purpose.py](https://codecov.io/gh/spacetx/starfish/pull/1371/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9pbWFnZXN0YWNrL3Rlc3QvZmFjdG9yaWVzL2FsbF9wdXJwb3NlLnB5) | `24.44% <0%> (-75.56%)` | :arrow_down: |
| [starfish/core/experiment/builder/\_\_init\_\_.py](https://codecov.io/gh/spacetx/starfish/pull/1371/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9leHBlcmltZW50L2J1aWxkZXIvX19pbml0X18ucHk=) | `30.35% <0%> (-69.65%)` | :arrow_down: |
| [starfish/core/spacetx\_format/test\_fuzz.py](https://codecov.io/gh/spacetx/starfish/pull/1371/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9zcGFjZXR4X2Zvcm1hdC90ZXN0X2Z1enoucHk=) | `31.03% <0%> (-68.97%)` | :arrow_down: |
| [starfish/core/intensity\_table/overlap.py](https://codecov.io/gh/spacetx/starfish/pull/1371/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9pbnRlbnNpdHlfdGFibGUvb3ZlcmxhcC5weQ==) | `31.37% <0%> (-64.71%)` | :arrow_down: |
| [.../spots/\_detect\_spots/local\_search\_blob\_detector.py](https://codecov.io/gh/spacetx/starfish/pull/1371/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9zcG90cy9fZGV0ZWN0X3Nwb3RzL2xvY2FsX3NlYXJjaF9ibG9iX2RldGVjdG9yLnB5) | `32.2% <0%> (-61.87%)` | :arrow_down: |
| ... and [65 more](https://codecov.io/gh/spacetx/starfish/pull/1371/diff?src=pr&el=tree-more) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/spacetx/starfish/pull/1371?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/1371?src=pr&el=footer). Last update [8af3e83...ebf9b14](https://codecov.io/gh/spacetx/starfish/pull/1371?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


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/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