From da247e65d8dbfaf1a5cfb8030b3e113048416b80 Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Fri, 1 Mar 2019 16:58:41 -0800 Subject: [PATCH 01/21] adding todo's and comment about looping --- ...peline_-_Human_Occipital_Cortex_-_1_FOV.py | 1 + starfish/experiment/experiment.py | 68 +++++++++---------- starfish/intensity_table/intensity_table.py | 9 ++- 3 files changed, 43 insertions(+), 35 deletions(-) diff --git a/notebooks/py/DARTFISH_Pipeline_-_Human_Occipital_Cortex_-_1_FOV.py b/notebooks/py/DARTFISH_Pipeline_-_Human_Occipital_Cortex_-_1_FOV.py index bd22392c6..472ff5e96 100644 --- a/notebooks/py/DARTFISH_Pipeline_-_Human_Occipital_Cortex_-_1_FOV.py +++ b/notebooks/py/DARTFISH_Pipeline_-_Human_Occipital_Cortex_-_1_FOV.py @@ -41,6 +41,7 @@ use_test_data = os.getenv("USE_TEST_DATA") is not None exp = data.DARTFISH(use_test_data=use_test_data) +# EXAMPLE SIMPLE LOOPING OVER ALIGNED IMAGE GROUPS: for stack in exp.fov().iterate_image_type(FieldOfView.PRIMARY_IMAGES) stack = exp.fov().get_image(FieldOfView.PRIMARY_IMAGES) # EPY: END code diff --git a/starfish/experiment/experiment.py b/starfish/experiment/experiment.py index e5e8b3b48..48914983d 100644 --- a/starfish/experiment/experiment.py +++ b/starfish/experiment/experiment.py @@ -1,11 +1,7 @@ -import copy import json -import pprint -from collections import OrderedDict from typing import ( Callable, Dict, - Iterator, List, MutableMapping, MutableSequence, @@ -39,8 +35,6 @@ class FieldOfView: type. The primary image is accessed using the name :py:attr:`starfish.experiment.experiment.FieldOFView.PRIMARY_IMAGES`. - Access a FOV through a experiment. experiement.fov() - Attributes ---------- name : str @@ -53,7 +47,7 @@ class FieldOfView: def __init__( self, name: str, - image_tilesets: MutableMapping[str, TileSet] + image_tilesets: Optional[MutableMapping[str, TileSet]] = None ) -> None: """ Fields of views can obtain their primary image from either an ImageStack or a TileSet (but @@ -68,9 +62,10 @@ def __init__( 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] = self.parse_coordinate_groups(tileset) - self._images = image_tilesets + if image_tilesets: + for name, tileset in image_tilesets.items(): + self.aligned_coordinate_groups[name] = self.parse_coordinate_groups(tileset) + self._images = image_tilesets def __repr__(self): images = '\n '.join( @@ -94,19 +89,18 @@ def parse_coordinate_groups(self, tileset: TileSet) -> List[CropParameters]: A list of CropParameters. Each entry describes the r/ch/z values of tiles that are aligned (have matching coordinates) """ - coord_groups: OrderedDict[tuple, CropParameters] = OrderedDict() - for tile in tileset.tiles(): - x_y_coords = ( + coord_groups: Dict[tuple, CropParameters] = dict() + for tile in tileset._tiles: + x_y_coords = tuple(i for i in [ 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 + tile.coordinates[Coordinates.Y][0], tile.coordinates[Coordinates.Y][1], + ]) 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]) + 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]) + crop_params.add_permitted_axes(Axes.ZPLANE, tile.indices[Axes.ZPLANE]) else: coord_groups[x_y_coords] = CropParameters( permitted_chs=[tile.indices[Axes.CH]], @@ -123,39 +117,45 @@ def name(self) -> str: def image_types(self) -> Set[str]: return set(self._images.keys()) - def show_aligned_image_groups(self) -> None: + def show_aligned_image_groups(self): """ Describe the aligned subgroups for each Tileset in this FOV ex. - {'nuclei': ' Group 0: ', - 'primary': ' Group 0: '} + {'primary': 'Group 0: ', + 'nuclei': 'Group 0: '} + + Means there are two tilesets in this FOV, (primary and nuclei) each tileset only + has one aligned subgroup. + + ex. + {'primary': 'Group 0: + Group 1: '} + + Means there is one tileset in this FOV, (primary), it has two different aligned + coordinate groups with 4 rounds in each - Means there are two tilesets in this FOV (primary and nuclei), and because all images have - the same (x, y) coordinates, each tileset has a single aligned subgroup. """ 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( + info = '\n '.join( f" Group {k}: " f" " + f"""({len(v._permitted_rounds) if v._permitted_rounds else 1, + len(v._permitted_chs) if v._permitted_chs else 1, + len(v._permitted_zplanes) if v._permitted_zplanes else 1, y_size, x_size})>""" for k, v in enumerate(groups) ) all_groups[name] = f'{info}' - pprint.pprint(all_groups) + return all_groups - def iterate_image_type(self, image_type: str) -> Iterator[ImageStack]: + def iterate_image_type(self, image_type: str): 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, aligned_group: int = 0, x_slice: Optional[Union[int, slice]] = None, y_slice: Optional[Union[int, slice]] = None, ) -> ImageStack: @@ -176,7 +176,7 @@ def get_image(self, item: str, aligned_group: int = 0, ------- The instantiated ImageStack """ - crop_params = copy.copy((self.aligned_coordinate_groups[item][aligned_group])) + crop_params = self.aligned_coordinate_groups[item][aligned_group] crop_params._x_slice = x_slice crop_params._y_slice = y_slice return ImageStack.from_tileset(self._images[item], crop_parameters=crop_params) diff --git a/starfish/intensity_table/intensity_table.py b/starfish/intensity_table/intensity_table.py index f9b849624..fef2f3fe8 100644 --- a/starfish/intensity_table/intensity_table.py +++ b/starfish/intensity_table/intensity_table.py @@ -1,6 +1,6 @@ from itertools import product from json import loads -from typing import Dict, Optional, Union +from typing import Dict, List, Optional, Union import numpy as np import pandas as pd @@ -388,6 +388,13 @@ def from_image_stack( return IntensityTable.from_spot_data(intensity_data, pixel_coordinates) + @staticmethod + def concatanate_intensity_tables(intensity_tables: List["IntensityTable"]): + # TODO VARY CONCAT LOGIC IF TILES OVERLAP + # This method is a starting point for handling tile overlap, right now + # it does a simple concat but people want other overlap logic implmented + return xr.concat(intensity_tables, dim=Features.AXIS) + def to_features_dataframe(self) -> pd.DataFrame: """Generates a dataframe of the underlying features multi-index. This is guaranteed to contain the features x, y, z, and radius. From d8122f596ae32b93f657ea850ec71570ab6a9f79 Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Fri, 1 Mar 2019 17:03:28 -0800 Subject: [PATCH 02/21] reverting experiement.py" --- starfish/experiment/experiment.py | 68 +++++++++++++++---------------- 1 file changed, 34 insertions(+), 34 deletions(-) diff --git a/starfish/experiment/experiment.py b/starfish/experiment/experiment.py index 48914983d..e5e8b3b48 100644 --- a/starfish/experiment/experiment.py +++ b/starfish/experiment/experiment.py @@ -1,7 +1,11 @@ +import copy import json +import pprint +from collections import OrderedDict from typing import ( Callable, Dict, + Iterator, List, MutableMapping, MutableSequence, @@ -35,6 +39,8 @@ class FieldOfView: type. The primary image is accessed using the name :py:attr:`starfish.experiment.experiment.FieldOFView.PRIMARY_IMAGES`. + Access a FOV through a experiment. experiement.fov() + Attributes ---------- name : str @@ -47,7 +53,7 @@ class FieldOfView: def __init__( self, name: str, - image_tilesets: Optional[MutableMapping[str, TileSet]] = None + image_tilesets: MutableMapping[str, TileSet] ) -> None: """ Fields of views can obtain their primary image from either an ImageStack or a TileSet (but @@ -62,10 +68,9 @@ def __init__( self._images: MutableMapping[str, TileSet] = dict() self._name = name self.aligned_coordinate_groups: Dict[str, List[CropParameters]] = dict() - if image_tilesets: - for name, tileset in image_tilesets.items(): - self.aligned_coordinate_groups[name] = self.parse_coordinate_groups(tileset) - self._images = image_tilesets + for name, tileset in image_tilesets.items(): + self.aligned_coordinate_groups[name] = self.parse_coordinate_groups(tileset) + self._images = image_tilesets def __repr__(self): images = '\n '.join( @@ -89,18 +94,19 @@ def parse_coordinate_groups(self, tileset: TileSet) -> List[CropParameters]: A list of CropParameters. Each entry describes the r/ch/z values of tiles that are aligned (have matching coordinates) """ - coord_groups: Dict[tuple, CropParameters] = dict() - for tile in tileset._tiles: - x_y_coords = tuple(i for i in [ + 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], - ]) + 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]) + 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]) + crop_params._add_permitted_axes(Axes.ZPLANE, tile.indices[Axes.ZPLANE]) else: coord_groups[x_y_coords] = CropParameters( permitted_chs=[tile.indices[Axes.CH]], @@ -117,45 +123,39 @@ def name(self) -> str: def image_types(self) -> Set[str]: return set(self._images.keys()) - def show_aligned_image_groups(self): + def show_aligned_image_groups(self) -> None: """ Describe the aligned subgroups for each Tileset in this FOV ex. - {'primary': 'Group 0: ', - 'nuclei': 'Group 0: '} - - Means there are two tilesets in this FOV, (primary and nuclei) each tileset only - has one aligned subgroup. - - ex. - {'primary': 'Group 0: - Group 1: '} - - Means there is one tileset in this FOV, (primary), it has two different aligned - coordinate groups with 4 rounds in each + {'nuclei': ' Group 0: ', + 'primary': ' Group 0: '} + Means there are two tilesets in this FOV (primary and nuclei), and because all images have + the same (x, y) coordinates, each tileset has a single aligned subgroup. """ 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( + info = '\n'.join( f" Group {k}: " f" """ + 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}' - return all_groups + pprint.pprint(all_groups) - def iterate_image_type(self, image_type: str): + def iterate_image_type(self, image_type: str) -> Iterator[ImageStack]: 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, aligned_group: int = 0, + def get_image(self, item: str, aligned_group: int = 0, x_slice: Optional[Union[int, slice]] = None, y_slice: Optional[Union[int, slice]] = None, ) -> ImageStack: @@ -176,7 +176,7 @@ def get_image(self, item, aligned_group: int = 0, ------- The instantiated ImageStack """ - crop_params = self.aligned_coordinate_groups[item][aligned_group] + crop_params = copy.copy((self.aligned_coordinate_groups[item][aligned_group])) crop_params._x_slice = x_slice crop_params._y_slice = y_slice return ImageStack.from_tileset(self._images[item], crop_parameters=crop_params) From 2fb6fa0aafa5ba9346876d993e02bc20ba680787 Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Mon, 4 Mar 2019 09:30:39 -0800 Subject: [PATCH 03/21] example looping on allen test --- ...peline_-_Human_Occipital_Cortex_-_1_FOV.py | 1 - .../full_pipelines/api/test_allen_smFISH.py | 62 +++++++++---------- 2 files changed, 31 insertions(+), 32 deletions(-) diff --git a/notebooks/py/DARTFISH_Pipeline_-_Human_Occipital_Cortex_-_1_FOV.py b/notebooks/py/DARTFISH_Pipeline_-_Human_Occipital_Cortex_-_1_FOV.py index 472ff5e96..bd22392c6 100644 --- a/notebooks/py/DARTFISH_Pipeline_-_Human_Occipital_Cortex_-_1_FOV.py +++ b/notebooks/py/DARTFISH_Pipeline_-_Human_Occipital_Cortex_-_1_FOV.py @@ -41,7 +41,6 @@ use_test_data = os.getenv("USE_TEST_DATA") is not None exp = data.DARTFISH(use_test_data=use_test_data) -# EXAMPLE SIMPLE LOOPING OVER ALIGNED IMAGE GROUPS: for stack in exp.fov().iterate_image_type(FieldOfView.PRIMARY_IMAGES) stack = exp.fov().get_image(FieldOfView.PRIMARY_IMAGES) # EPY: END code diff --git a/starfish/test/full_pipelines/api/test_allen_smFISH.py b/starfish/test/full_pipelines/api/test_allen_smFISH.py index b888bd222..327f44b01 100644 --- a/starfish/test/full_pipelines/api/test_allen_smFISH.py +++ b/starfish/test/full_pipelines/api/test_allen_smFISH.py @@ -1,15 +1,13 @@ import numpy as np -import pytest import starfish.data -from starfish import FieldOfView +from starfish import FieldOfView, IntensityTable from starfish.image._filter.bandpass import Bandpass from starfish.image._filter.clip import Clip from starfish.image._filter.gaussian_low_pass import GaussianLowPass from starfish.spots._detector.trackpy_local_max_peak_finder import TrackpyLocalMaxPeakFinder -@pytest.mark.skip('issues with checksums prevent this data from working properly') def test_allen_smFISH_cropped_data(): # set random seed to errors provoked by optimization functions @@ -17,31 +15,33 @@ def test_allen_smFISH_cropped_data(): # load the experiment experiment = starfish.data.allen_smFISH(use_test_data=True) - - primary_image = experiment.fov().get_image(FieldOfView.PRIMARY_IMAGES) - - clip = Clip(p_min=10, p_max=100) - clipped_image = clip.run(primary_image, in_place=False) - - bandpass = Bandpass(lshort=0.5, llong=7, threshold=None, truncate=4) - bandpassed_image = bandpass.run(clipped_image, in_place=False) - - clip = Clip(p_min=10, p_max=100, is_volume=False) - clipped_bandpassed_image = clip.run(bandpassed_image, in_place=False) - - sigma = (1, 0, 0) # filter only in z, do nothing in x, y - glp = GaussianLowPass(sigma=sigma, is_volume=True) - z_filtered_image = glp.run(clipped_bandpassed_image, in_place=False) - - lmpf = TrackpyLocalMaxPeakFinder( - spot_diameter=3, - min_mass=300, - max_size=3, - separation=5, - noise_size=0.65, - preprocess=False, - percentile=10, - verbose=True, - is_volume=True, - ) - intensities = lmpf.run(z_filtered_image) # noqa + all_intensities = list() + for primary_image in experiment.fov().iterate_image_type(FieldOfView.PRIMARY_IMAGES): + + clip = Clip(p_min=10, p_max=100) + clipped_image = clip.run(primary_image, in_place=False) + + bandpass = Bandpass(lshort=0.5, llong=7, truncate=4) + bandpassed_image = bandpass.run(clipped_image, in_place=False) + + clip = Clip(p_min=10, p_max=100, is_volume=False) + clipped_bandpassed_image = clip.run(bandpassed_image, in_place=False) + + sigma = (1, 0, 0) # filter only in z, do nothing in x, y + glp = GaussianLowPass(sigma=sigma, is_volume=True) + z_filtered_image = glp.run(clipped_bandpassed_image, in_place=False) + + lmpf = TrackpyLocalMaxPeakFinder( + spot_diameter=3, + min_mass=300, + max_size=3, + separation=5, + noise_size=0.65, + preprocess=False, + percentile=10, + verbose=True, + is_volume=True, + ) + intensities = lmpf.run(z_filtered_image) # noqa + all_intensities.append(intensities) + IntensityTable.concatanate_intensity_tables(all_intensities) From cb17b5f2560e67cad27d946f08a41f5cfc07ec0c Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Mon, 4 Mar 2019 10:28:25 -0800 Subject: [PATCH 04/21] adding travis wait --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index f6e0a2a4f..0fcc2c6da 100644 --- a/.travis.yml +++ b/.travis.yml @@ -12,7 +12,7 @@ stages: jobs: include: - name: Syntax and Fast Tests - script: make install-dev fast + script: travis_wait make install-dev fast after_success: - bash <(curl -s https://codecov.io/bash) - name: Slow Tests From 7de5bd633adeb035ae6cdeb155ad9c33155ad757 Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Mon, 4 Mar 2019 11:43:24 -0800 Subject: [PATCH 05/21] upping travis wait time --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 0fcc2c6da..14bda885b 100644 --- a/.travis.yml +++ b/.travis.yml @@ -12,7 +12,7 @@ stages: jobs: include: - name: Syntax and Fast Tests - script: travis_wait make install-dev fast + script: travis_wait 30 make install-dev fast after_success: - bash <(curl -s https://codecov.io/bash) - name: Slow Tests From 4db58d998ab2fb636ebb4eaf5e74cbf9c7db855e Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Tue, 5 Mar 2019 11:58:45 -0800 Subject: [PATCH 06/21] trying to free up mem in allen test --- starfish/test/full_pipelines/api/test_allen_smFISH.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/starfish/test/full_pipelines/api/test_allen_smFISH.py b/starfish/test/full_pipelines/api/test_allen_smFISH.py index 327f44b01..9c3c03d12 100644 --- a/starfish/test/full_pipelines/api/test_allen_smFISH.py +++ b/starfish/test/full_pipelines/api/test_allen_smFISH.py @@ -1,3 +1,5 @@ +import gc + import numpy as np import starfish.data @@ -44,4 +46,5 @@ def test_allen_smFISH_cropped_data(): ) intensities = lmpf.run(z_filtered_image) # noqa all_intensities.append(intensities) + gc.collect() IntensityTable.concatanate_intensity_tables(all_intensities) From 5a90da732e0562d8b4883d0d7bf5b66f9801732b Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Tue, 5 Mar 2019 12:02:26 -0800 Subject: [PATCH 07/21] reverting mem thing --- starfish/test/full_pipelines/api/test_allen_smFISH.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/starfish/test/full_pipelines/api/test_allen_smFISH.py b/starfish/test/full_pipelines/api/test_allen_smFISH.py index 9c3c03d12..327f44b01 100644 --- a/starfish/test/full_pipelines/api/test_allen_smFISH.py +++ b/starfish/test/full_pipelines/api/test_allen_smFISH.py @@ -1,5 +1,3 @@ -import gc - import numpy as np import starfish.data @@ -46,5 +44,4 @@ def test_allen_smFISH_cropped_data(): ) intensities = lmpf.run(z_filtered_image) # noqa all_intensities.append(intensities) - gc.collect() IntensityTable.concatanate_intensity_tables(all_intensities) From a9fb8a96a5cd566eb92ac5e2be5e4ddce291dc9c Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Tue, 5 Mar 2019 12:58:18 -0800 Subject: [PATCH 08/21] moving ecample to 3dsmFish --- notebooks/py/3d_smFISH.py | 36 ++++++----- .../full_pipelines/api/test_allen_smFISH.py | 62 +++++++++---------- 2 files changed, 51 insertions(+), 47 deletions(-) diff --git a/notebooks/py/3d_smFISH.py b/notebooks/py/3d_smFISH.py index 402490e31..10970bbfc 100644 --- a/notebooks/py/3d_smFISH.py +++ b/notebooks/py/3d_smFISH.py @@ -30,7 +30,8 @@ import starfish import starfish.display -from starfish import data, FieldOfView +from starfish import data, FieldOfView, IntensityTable + # EPY: END code # EPY: START markdown @@ -99,22 +100,25 @@ def processing_pipeline( print("Loading images...") primary_image = experiment[fov_name].get_image(FieldOfView.PRIMARY_IMAGES) + all_intensities = list() codebook = experiment.codebook - - print("Filtering images...") - filter_kwargs = dict( - in_place=True, - verbose=True, - n_processes=n_processes - ) - clip1.run(primary_image, **filter_kwargs) - bandpass.run(primary_image, **filter_kwargs) - glp.run(primary_image, **filter_kwargs) - clip2.run(primary_image, **filter_kwargs) - - print("Calling spots...") - spot_attributes = tlmpf.run(primary_image) - + for primary_image in experiment[fov_name].iterate_image_type(FieldOfView.PRIMARY_IMAGES): + print("Filtering images...") + filter_kwargs = dict( + in_place=True, + verbose=True, + n_processes=n_processes + ) + clip1.run(primary_image, **filter_kwargs) + bandpass.run(primary_image, **filter_kwargs) + glp.run(primary_image, **filter_kwargs) + clip2.run(primary_image, **filter_kwargs) + + print("Calling spots...") + spot_attributes = tlmpf.run(primary_image) + all_intensities.append(spot_attributes) + + spot_attributes = IntensityTable.concatanate_intensity_tables(all_intensities) print("Decoding spots...") decoded = codebook.decode_per_round_max(spot_attributes) decoded = decoded[decoded["total_intensity"]>.025] diff --git a/starfish/test/full_pipelines/api/test_allen_smFISH.py b/starfish/test/full_pipelines/api/test_allen_smFISH.py index 327f44b01..b888bd222 100644 --- a/starfish/test/full_pipelines/api/test_allen_smFISH.py +++ b/starfish/test/full_pipelines/api/test_allen_smFISH.py @@ -1,13 +1,15 @@ import numpy as np +import pytest import starfish.data -from starfish import FieldOfView, IntensityTable +from starfish import FieldOfView from starfish.image._filter.bandpass import Bandpass from starfish.image._filter.clip import Clip from starfish.image._filter.gaussian_low_pass import GaussianLowPass from starfish.spots._detector.trackpy_local_max_peak_finder import TrackpyLocalMaxPeakFinder +@pytest.mark.skip('issues with checksums prevent this data from working properly') def test_allen_smFISH_cropped_data(): # set random seed to errors provoked by optimization functions @@ -15,33 +17,31 @@ def test_allen_smFISH_cropped_data(): # load the experiment experiment = starfish.data.allen_smFISH(use_test_data=True) - all_intensities = list() - for primary_image in experiment.fov().iterate_image_type(FieldOfView.PRIMARY_IMAGES): - - clip = Clip(p_min=10, p_max=100) - clipped_image = clip.run(primary_image, in_place=False) - - bandpass = Bandpass(lshort=0.5, llong=7, truncate=4) - bandpassed_image = bandpass.run(clipped_image, in_place=False) - - clip = Clip(p_min=10, p_max=100, is_volume=False) - clipped_bandpassed_image = clip.run(bandpassed_image, in_place=False) - - sigma = (1, 0, 0) # filter only in z, do nothing in x, y - glp = GaussianLowPass(sigma=sigma, is_volume=True) - z_filtered_image = glp.run(clipped_bandpassed_image, in_place=False) - - lmpf = TrackpyLocalMaxPeakFinder( - spot_diameter=3, - min_mass=300, - max_size=3, - separation=5, - noise_size=0.65, - preprocess=False, - percentile=10, - verbose=True, - is_volume=True, - ) - intensities = lmpf.run(z_filtered_image) # noqa - all_intensities.append(intensities) - IntensityTable.concatanate_intensity_tables(all_intensities) + + primary_image = experiment.fov().get_image(FieldOfView.PRIMARY_IMAGES) + + clip = Clip(p_min=10, p_max=100) + clipped_image = clip.run(primary_image, in_place=False) + + bandpass = Bandpass(lshort=0.5, llong=7, threshold=None, truncate=4) + bandpassed_image = bandpass.run(clipped_image, in_place=False) + + clip = Clip(p_min=10, p_max=100, is_volume=False) + clipped_bandpassed_image = clip.run(bandpassed_image, in_place=False) + + sigma = (1, 0, 0) # filter only in z, do nothing in x, y + glp = GaussianLowPass(sigma=sigma, is_volume=True) + z_filtered_image = glp.run(clipped_bandpassed_image, in_place=False) + + lmpf = TrackpyLocalMaxPeakFinder( + spot_diameter=3, + min_mass=300, + max_size=3, + separation=5, + noise_size=0.65, + preprocess=False, + percentile=10, + verbose=True, + is_volume=True, + ) + intensities = lmpf.run(z_filtered_image) # noqa From 9816592dfa8eff5b3030522d2febd9ab3637e205 Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Tue, 5 Mar 2019 13:00:01 -0800 Subject: [PATCH 09/21] reverting tavis wait --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 14bda885b..f6e0a2a4f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -12,7 +12,7 @@ stages: jobs: include: - name: Syntax and Fast Tests - script: travis_wait 30 make install-dev fast + script: make install-dev fast after_success: - bash <(curl -s https://codecov.io/bash) - name: Slow Tests From 84b00f61f212984c449aae130a7d8a614239ea57 Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Tue, 5 Mar 2019 13:38:01 -0800 Subject: [PATCH 10/21] regenerated notebooks --- notebooks/3d_smFISH.ipynb | 36 ++++++++++++++++++-------------- notebooks/assay_comparison.ipynb | 2 +- 2 files changed, 21 insertions(+), 17 deletions(-) diff --git a/notebooks/3d_smFISH.ipynb b/notebooks/3d_smFISH.ipynb index 03f41fa61..e78168885 100644 --- a/notebooks/3d_smFISH.ipynb +++ b/notebooks/3d_smFISH.ipynb @@ -41,7 +41,8 @@ "\n", "import starfish\n", "import starfish.display\n", - "from starfish import data, FieldOfView" + "from starfish import data, FieldOfView, IntensityTable\n", + "" ] }, { @@ -129,22 +130,25 @@ "\n", " print(\"Loading images...\")\n", " primary_image = experiment[fov_name].get_image(FieldOfView.PRIMARY_IMAGES)\n", + " all_intensities = list()\n", " codebook = experiment.codebook\n", - "\n", - " print(\"Filtering images...\")\n", - " filter_kwargs = dict(\n", - " in_place=True,\n", - " verbose=True,\n", - " n_processes=n_processes\n", - " )\n", - " clip1.run(primary_image, **filter_kwargs)\n", - " bandpass.run(primary_image, **filter_kwargs)\n", - " glp.run(primary_image, **filter_kwargs)\n", - " clip2.run(primary_image, **filter_kwargs)\n", - "\n", - " print(\"Calling spots...\")\n", - " spot_attributes = tlmpf.run(primary_image)\n", - "\n", + " for primary_image in experiment[fov_name].iterate_image_type(FieldOfView.PRIMARY_IMAGES):\n", + " print(\"Filtering images...\")\n", + " filter_kwargs = dict(\n", + " in_place=True,\n", + " verbose=True,\n", + " n_processes=n_processes\n", + " )\n", + " clip1.run(primary_image, **filter_kwargs)\n", + " bandpass.run(primary_image, **filter_kwargs)\n", + " glp.run(primary_image, **filter_kwargs)\n", + " clip2.run(primary_image, **filter_kwargs)\n", + "\n", + " print(\"Calling spots...\")\n", + " spot_attributes = tlmpf.run(primary_image)\n", + " all_intensities.append(spot_attributes)\n", + "\n", + " spot_attributes = IntensityTable.concatanate_intensity_tables(all_intensities)\n", " print(\"Decoding spots...\")\n", " decoded = codebook.decode_per_round_max(spot_attributes)\n", " decoded = decoded[decoded[\"total_intensity\"]>.025]\n", diff --git a/notebooks/assay_comparison.ipynb b/notebooks/assay_comparison.ipynb index 2b29064f8..5fe797bc4 100644 --- a/notebooks/assay_comparison.ipynb +++ b/notebooks/assay_comparison.ipynb @@ -482,4 +482,4 @@ }, "nbformat": 4, "nbformat_minor": 2 -} +} \ No newline at end of file From fb1764d63ef6be5ed3ff8a507c8a0d19ad214fa4 Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Thu, 7 Mar 2019 09:32:15 -0800 Subject: [PATCH 11/21] :changes --- starfish/intensity_table/intensity_table.py | 43 +++++++++++++++++++-- starfish/types/_constants.py | 4 ++ 2 files changed, 43 insertions(+), 4 deletions(-) diff --git a/starfish/intensity_table/intensity_table.py b/starfish/intensity_table/intensity_table.py index fef2f3fe8..d990e8c15 100644 --- a/starfish/intensity_table/intensity_table.py +++ b/starfish/intensity_table/intensity_table.py @@ -10,6 +10,7 @@ from starfish.expression_matrix.expression_matrix import ExpressionMatrix from starfish.image._filter.util import preserve_float_range from starfish.types import Axes, Features, LOG, SpotAttributes, STARFISH_EXTRAS_KEY +from starfish.types._constants import OverlapStrategy class IntensityTable(xr.DataArray): @@ -389,10 +390,44 @@ def from_image_stack( return IntensityTable.from_spot_data(intensity_data, pixel_coordinates) @staticmethod - def concatanate_intensity_tables(intensity_tables: List["IntensityTable"]): - # TODO VARY CONCAT LOGIC IF TILES OVERLAP - # This method is a starting point for handling tile overlap, right now - # it does a simple concat but people want other overlap logic implmented + def take_max(intensity_table1: "IntensityTable", intensityTable2: "IntensityTable") -> None: + """Find the overlapping sections between IntensityTables + and eliminate spots from sections with less spots. + """ + + # Find intersection of two arrays + i1_min_x, i1_max_x, i1_min_y, i1_max_y = min(intensity_table1['xc']).data, \ + max(intensity_table1['xc']).data, \ + min(intensity_table1['yc']).data,\ + max(intensity_table1['yc']).data + + i2_min_x, i2_max_x, i2_min_y, i2_max_y = min(intensity_table1['xc']).data, \ + max(intensity_table1['xc']).data, \ + min(intensity_table1['yc']).data, \ + max(intensity_table1['yc']).data + + # select from i1 where coord values in 12min/max + + # select from 12 where coord values in i1 min/mx + + + + # Count spots in intersection for both + # Whichever has the least spots, convert values to nan + + + return None + + OVERLAP_STRATEGY_MAP = { + OverlapStrategy.TAKE_MAX: take_max + } + + @staticmethod + def concatanate_intensity_tables(intensity_tables: List["IntensityTable"], + overlap_strategy: OverlapStrategy = None): + if overlap_strategy: + overlap_method = IntensityTable.OVERLAP_STRATEGY_MAP[overlap_strategy] + overlap_method(intensity_tables) return xr.concat(intensity_tables, dim=Features.AXIS) def to_features_dataframe(self) -> pd.DataFrame: diff --git a/starfish/types/_constants.py b/starfish/types/_constants.py index aea2e4c6b..a3f5fecaa 100644 --- a/starfish/types/_constants.py +++ b/starfish/types/_constants.py @@ -76,3 +76,7 @@ class Features: CELL_ID = 'cell_id' SPOT_ID = 'spot_id' INTENSITY = 'intensity' + + +class OverlapStrategy: + TAKE_MAX = 'take_max' From 29ff4b8668272e7db0e99f50c83cf3cc9ec034f8 Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Fri, 22 Mar 2019 15:57:55 -0700 Subject: [PATCH 12/21] changes. --- starfish/intensity_table/intensity_table.py | 48 ++++++++++++------- .../test_intensity_table_concat.py | 33 +++++++++++++ 2 files changed, 65 insertions(+), 16 deletions(-) create mode 100644 starfish/test/intensity_table/test_intensity_table_concat.py diff --git a/starfish/intensity_table/intensity_table.py b/starfish/intensity_table/intensity_table.py index a94a058f0..9092aeea9 100644 --- a/starfish/intensity_table/intensity_table.py +++ b/starfish/intensity_table/intensity_table.py @@ -390,34 +390,50 @@ def from_image_stack( return IntensityTable.from_spot_data(intensity_data, pixel_coordinates) + + @staticmethod - def take_max(intensity_table1: "IntensityTable", intensityTable2: "IntensityTable") -> None: + def take_max(i1: "IntensityTable", i2: "IntensityTable") -> None: """Find the overlapping sections between IntensityTables - and eliminate spots from sections with less spots. + and eliminate spots from whichever section has less spots. """ - # Find intersection of two arrays - i1_min_x, i1_max_x, i1_min_y, i1_max_y = min(intensity_table1['xc']).data, \ - max(intensity_table1['xc']).data, \ - min(intensity_table1['yc']).data,\ - max(intensity_table1['yc']).data + # Get bounding physical coords + i1_min_x, i1_max_x, i1_min_y, i1_max_y = min(i1['xc']).data, \ + max(i1['xc']).data, \ + min(i1['yc']).data,\ + max(i1['yc']).data - i2_min_x, i2_max_x, i2_min_y, i2_max_y = min(intensity_table1['xc']).data, \ - max(intensity_table1['xc']).data, \ - min(intensity_table1['yc']).data, \ - max(intensity_table1['yc']).data + i2_min_x, i2_max_x, i2_min_y, i2_max_y = min(i2['xc']).data, \ + max(i2['xc']).data, \ + min(i2['yc']).data, \ + max(i2['yc']).data - # select from i1 where coord values in 12min/max - # select from 12 where coord values in i1 min/mx + overlap_xmin = max(i1_min_x, i2_min_x) + overlap_xmax = min(i1_max_x, i2_max_y) + overlap_ymin = max(i1_min_y, i2_min_y) + overlap_ymax = min(i1_max_y, i2_max_y) + intersect1 = i1.where((i1.xc > overlap_xmin) & + (i1.xc < overlap_xmax) & + (i1.yc > overlap_ymin) & + (i1.yc < overlap_ymax), drop=True) + + # select from 12 where coord values in i1 min/mx + intersect2 = i1.where((i1.xc > overlap_xmin) & + (i1.xc < overlap_xmax) & + (i1.yc > overlap_ymin) & + (i1.yc < overlap_ymax), drop=True) - # Count spots in intersection for both - # Whichever has the least spots, convert values to nan + if intersect1.sizes[Features.AXIS] > intersect2.sizes[Features.AXIS]: + # zero out spots in intersect1 + else: + # zero out spots in interset2 - return None + return intersect1 OVERLAP_STRATEGY_MAP = { OverlapStrategy.TAKE_MAX: take_max diff --git a/starfish/test/intensity_table/test_intensity_table_concat.py b/starfish/test/intensity_table/test_intensity_table_concat.py new file mode 100644 index 000000000..baaef28a8 --- /dev/null +++ b/starfish/test/intensity_table/test_intensity_table_concat.py @@ -0,0 +1,33 @@ +import numpy as np +import xarray as xr + +from starfish import IntensityTable +from starfish.test import test_utils +from starfish.types import Coordinates + + +def test_overlap(): + codebook = test_utils.codebook_array_factory() + it1 = IntensityTable.synthetic_intensities( + codebook, + num_z=1, + height=50, + width=50, + n_spots=10 + ) + # intensity table 1 has 10 spots, xmin = 0, ymin = 0, xmax = 2, ymax = 1 + it1[Coordinates.X.value] = xr.DataArray(np.linspace(0, 2, 10), dims='features') + it1[Coordinates.Y.value] = xr.DataArray(np.linspace(0, 1, 10), dims='features') + + it2 = IntensityTable.synthetic_intensities( + codebook, + num_z=1, + height=50, + width=50, + n_spots=12 + ) + # intensity table 2 has 12 spots, xmin = 1, ymin = 1, xmax = 3, ymax = 3 + it2[Coordinates.X.value] = xr.DataArray(np.linspace(1, 3, 12), dims='features') + it2[Coordinates.Y.value] = xr.DataArray(np.linspace(1, 3, 12), dims='features') + + IntensityTable.take_max(it1, it2) From b20c92c21858ecc310ad8d24a899c14f78cd0147 Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Wed, 27 Mar 2019 15:34:44 -0700 Subject: [PATCH 13/21] good framework to process lists of tables --- starfish/intensity_table/intensity_table.py | 72 +++++------- .../test_intensity_table_concat.py | 20 ++-- starfish/types/_constants.py | 1 + starfish/util/geometry.py | 103 ++++++++++++++++++ 4 files changed, 143 insertions(+), 53 deletions(-) create mode 100644 starfish/util/geometry.py diff --git a/starfish/intensity_table/intensity_table.py b/starfish/intensity_table/intensity_table.py index 9092aeea9..eca332e41 100644 --- a/starfish/intensity_table/intensity_table.py +++ b/starfish/intensity_table/intensity_table.py @@ -12,6 +12,12 @@ from starfish.types._constants import OverlapStrategy from starfish.types import Axes, DecodedSpots, Features, LOG, SpotAttributes, STARFISH_EXTRAS_KEY from starfish.util.dtype import preserve_float_range +from starfish.util.geometry import( + find_overlaps_of_xarrays, + OVERLAP_STRATEGY_MAP, + sel_area_of_xarray, + take_max +) class IntensityTable(xr.DataArray): @@ -393,59 +399,33 @@ def from_image_stack( @staticmethod - def take_max(i1: "IntensityTable", i2: "IntensityTable") -> None: + def process_overlaps(its: List["IntensityTable"], + overlap_strategy: OverlapStrategy + ) -> List["IntensityTable"]: """Find the overlapping sections between IntensityTables - and eliminate spots from whichever section has less spots. + and process them according to the give overlap strategy """ + overlap_pairs = find_overlaps_of_xarrays(its) + for idx1, idx2, intersection_rect in overlap_pairs: + intersect1 = sel_area_of_xarray(its[idx1], intersection_rect) + intersect2 = sel_area_of_xarray(its[idx2], intersection_rect) + # compare + overlap_mathod = OVERLAP_STRATEGY_MAP[overlap_strategy] + it1, it2 = overlap_mathod(intersection_rect, + its[idx1], intersect1, + its[idx2], intersect2) + its[idx1] = it1 + its[idx2] = it2 + return its - # Get bounding physical coords - i1_min_x, i1_max_x, i1_min_y, i1_max_y = min(i1['xc']).data, \ - max(i1['xc']).data, \ - min(i1['yc']).data,\ - max(i1['yc']).data - - i2_min_x, i2_max_x, i2_min_y, i2_max_y = min(i2['xc']).data, \ - max(i2['xc']).data, \ - min(i2['yc']).data, \ - max(i2['yc']).data - - - overlap_xmin = max(i1_min_x, i2_min_x) - overlap_xmax = min(i1_max_x, i2_max_y) - overlap_ymin = max(i1_min_y, i2_min_y) - overlap_ymax = min(i1_max_y, i2_max_y) - - intersect1 = i1.where((i1.xc > overlap_xmin) & - (i1.xc < overlap_xmax) & - (i1.yc > overlap_ymin) & - (i1.yc < overlap_ymax), drop=True) - - # select from 12 where coord values in i1 min/mx - intersect2 = i1.where((i1.xc > overlap_xmin) & - (i1.xc < overlap_xmax) & - (i1.yc > overlap_ymin) & - (i1.yc < overlap_ymax), drop=True) - - - - if intersect1.sizes[Features.AXIS] > intersect2.sizes[Features.AXIS]: - # zero out spots in intersect1 - else: - # zero out spots in interset2 - - return intersect1 - - OVERLAP_STRATEGY_MAP = { - OverlapStrategy.TAKE_MAX: take_max - } @staticmethod def concatanate_intensity_tables(intensity_tables: List["IntensityTable"], + process_overlaps: bool = False, overlap_strategy: OverlapStrategy = None): - if overlap_strategy: - overlap_method = IntensityTable.OVERLAP_STRATEGY_MAP[overlap_strategy] - overlap_method(intensity_tables) - + if process_overlaps: + intensity_tables = IntensityTable.process_overlaps(intensity_tables, + overlap_strategy) return xr.concat(intensity_tables, dim=Features.AXIS) def to_features_dataframe(self) -> pd.DataFrame: diff --git a/starfish/test/intensity_table/test_intensity_table_concat.py b/starfish/test/intensity_table/test_intensity_table_concat.py index baaef28a8..73ca8a58d 100644 --- a/starfish/test/intensity_table/test_intensity_table_concat.py +++ b/starfish/test/intensity_table/test_intensity_table_concat.py @@ -3,7 +3,8 @@ from starfish import IntensityTable from starfish.test import test_utils -from starfish.types import Coordinates +from starfish.types import Coordinates, Features +from starfish.types._constants import OverlapStrategy def test_overlap(): @@ -17,17 +18,22 @@ def test_overlap(): ) # intensity table 1 has 10 spots, xmin = 0, ymin = 0, xmax = 2, ymax = 1 it1[Coordinates.X.value] = xr.DataArray(np.linspace(0, 2, 10), dims='features') - it1[Coordinates.Y.value] = xr.DataArray(np.linspace(0, 1, 10), dims='features') + it1[Coordinates.Y.value] = xr.DataArray(np.linspace(0, 2, 10), dims='features') it2 = IntensityTable.synthetic_intensities( codebook, num_z=1, height=50, width=50, - n_spots=12 + n_spots=20 ) - # intensity table 2 has 12 spots, xmin = 1, ymin = 1, xmax = 3, ymax = 3 - it2[Coordinates.X.value] = xr.DataArray(np.linspace(1, 3, 12), dims='features') - it2[Coordinates.Y.value] = xr.DataArray(np.linspace(1, 3, 12), dims='features') + # intensity table 2 has 20 spots, xmin = 1, ymin = 1, xmax = 3, ymax = 3 + it2[Coordinates.X.value] = xr.DataArray(np.linspace(1, 3, 20), dims='features') + it2[Coordinates.Y.value] = xr.DataArray(np.linspace(1, 3, 20), dims='features') + + concatenated = IntensityTable.concatanate_intensity_tables([it1, it2], + process_overlaps=True, + overlap_strategy=OverlapStrategy.TAKE_MAX) + + assert concatenated.sizes[Features.AXIS] == 26 - IntensityTable.take_max(it1, it2) diff --git a/starfish/types/_constants.py b/starfish/types/_constants.py index 3f6015f1d..4315d2f8b 100644 --- a/starfish/types/_constants.py +++ b/starfish/types/_constants.py @@ -82,6 +82,7 @@ class OverlapStrategy: TAKE_MAX = 'take_max' + class Clip(AugmentedEnum): """ contains clipping options that determine how out-of-bounds values produced by filters are diff --git a/starfish/util/geometry.py b/starfish/util/geometry.py new file mode 100644 index 000000000..2c455f428 --- /dev/null +++ b/starfish/util/geometry.py @@ -0,0 +1,103 @@ +import itertools +from typing import List, Union, Tuple +import xarray as xr + +from starfish.types import Features +from starfish.types._constants import OverlapStrategy + + +class Rect: + def __init__(self, min_x, max_x, min_y, max_y): + """ + Small class that defines a rectangle by its bottom left + and top right coordinates. + """ + self.min_x = min_x + self.max_x = max_x + self.min_y = min_y + self.max_y = max_y + + @staticmethod + def no_overlap(rect1: "Rect", rect2: "Rect") -> bool: + """Return True if two rectangles do not overlap""" + return (rect1.max_x < rect2.min_x) | \ + (rect1.min_x > rect2.max_x) | \ + (rect1.max_y < rect2.min_y) | \ + (rect1.min_y > rect2.max_y) + + @staticmethod + def find_intersection(rect1: "Rect", rect2: "Rect") -> Union[None, "Rect"]: + """ + Find the intersection area of two intensity tables and return as Rect object. + If no overlap return none. + """ + if Rect.no_overlap(rect1, rect2): + return None + return Rect(max(rect1.min_x, rect2.min_x), + min(rect1.max_x, rect2.max_x), + max(rect1.min_y, rect2.min_y), + min(rect1.max_y, rect2.max_y)) + + +def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] + ) -> List[Tuple[int, int, "Rect"]]: + """ + Takes a list of xarrays and returns a list of xarray pairs that overlap physically + and their areas of overlap defined by a Rectangle. + """ + all_overlaps: List[Tuple[int, int, "Rect"]] = list() + for idx1, idx2 in itertools.combinations(range(len(xarrays)), 2): + xr1 = xarrays[idx1] + xr2 = xarrays[idx2] + rect1 = Rect(min(xr1['xc']).data, + max(xr1['xc']).data, + min(xr1['yc']).data, + max(xr1['yc']).data) + + rect2 = Rect(min(xr2['xc']).data, + max(xr2['xc']).data, + min(xr2['yc']).data, + max(xr2['yc']).data) + intersection = Rect.find_intersection(rect1, rect2) + if intersection: + all_overlaps.append((idx1, idx2, intersection)) + return all_overlaps + + +def remove_area_of_xarray(it: xr.DataArray, area: Rect) -> xr.DataArray: + """Remove section of xarray within a defined area""" + return it.where((it.xc <= area.min_x) | + (it.xc >= area.max_x) | + (it.yc <= area.min_y) | + (it.yc >= area.max_y), drop=True) + + +def sel_area_of_xarray(it: xr.DataArray, area: Rect) -> xr.DataArray: + """Select on xarray within a defined area""" + return it.where((it.xc > area.min_x) & + (it.xc < area.max_x) & + (it.yc > area.min_y) & + (it.yc < area.max_y), drop=True) + + +def take_max(intersection_rect: Rect, + it1: xr.DataArray, + intersect1: xr.DataArray, + it2: xr.DataArray, + intersect2: xr.DataArray): + # # compare to see which section has more spots + if intersect1.sizes[Features.AXIS] >= intersect2.sizes[Features.AXIS]: + # I1 has more spots remove intesection2 from it2 + it2 = remove_area_of_xarray(it2, intersection_rect) + else: + it1 = remove_area_of_xarray(it1, intersection_rect) + return it1, it2 + + +OVERLAP_STRATEGY_MAP = { + OverlapStrategy.TAKE_MAX: take_max +} + + + + From c52310c61d14fdea6aa27b3bffc000cb191dea78 Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Wed, 27 Mar 2019 15:51:28 -0700 Subject: [PATCH 14/21] cleanned up --- starfish/intensity_table/intensity_table.py | 17 +--- starfish/util/geometry.py | 103 ------------------- starfish/util/overlap_utils.py | 104 ++++++++++++++++++++ 3 files changed, 108 insertions(+), 116 deletions(-) delete mode 100644 starfish/util/geometry.py create mode 100644 starfish/util/overlap_utils.py diff --git a/starfish/intensity_table/intensity_table.py b/starfish/intensity_table/intensity_table.py index eca332e41..a5ba23de6 100644 --- a/starfish/intensity_table/intensity_table.py +++ b/starfish/intensity_table/intensity_table.py @@ -12,11 +12,9 @@ from starfish.types._constants import OverlapStrategy from starfish.types import Axes, DecodedSpots, Features, LOG, SpotAttributes, STARFISH_EXTRAS_KEY from starfish.util.dtype import preserve_float_range -from starfish.util.geometry import( +from starfish.util.overlap_utils import( find_overlaps_of_xarrays, OVERLAP_STRATEGY_MAP, - sel_area_of_xarray, - take_max ) @@ -396,8 +394,6 @@ def from_image_stack( return IntensityTable.from_spot_data(intensity_data, pixel_coordinates) - - @staticmethod def process_overlaps(its: List["IntensityTable"], overlap_strategy: OverlapStrategy @@ -407,18 +403,13 @@ def process_overlaps(its: List["IntensityTable"], """ overlap_pairs = find_overlaps_of_xarrays(its) for idx1, idx2, intersection_rect in overlap_pairs: - intersect1 = sel_area_of_xarray(its[idx1], intersection_rect) - intersect2 = sel_area_of_xarray(its[idx2], intersection_rect) - # compare - overlap_mathod = OVERLAP_STRATEGY_MAP[overlap_strategy] - it1, it2 = overlap_mathod(intersection_rect, - its[idx1], intersect1, - its[idx2], intersect2) + overlap_method = OVERLAP_STRATEGY_MAP[overlap_strategy] + # modify intensity tables based on overlap strategy + it1, it2 = overlap_method(intersection_rect, its[idx1], its[idx2]) its[idx1] = it1 its[idx2] = it2 return its - @staticmethod def concatanate_intensity_tables(intensity_tables: List["IntensityTable"], process_overlaps: bool = False, diff --git a/starfish/util/geometry.py b/starfish/util/geometry.py deleted file mode 100644 index 2c455f428..000000000 --- a/starfish/util/geometry.py +++ /dev/null @@ -1,103 +0,0 @@ -import itertools -from typing import List, Union, Tuple -import xarray as xr - -from starfish.types import Features -from starfish.types._constants import OverlapStrategy - - -class Rect: - def __init__(self, min_x, max_x, min_y, max_y): - """ - Small class that defines a rectangle by its bottom left - and top right coordinates. - """ - self.min_x = min_x - self.max_x = max_x - self.min_y = min_y - self.max_y = max_y - - @staticmethod - def no_overlap(rect1: "Rect", rect2: "Rect") -> bool: - """Return True if two rectangles do not overlap""" - return (rect1.max_x < rect2.min_x) | \ - (rect1.min_x > rect2.max_x) | \ - (rect1.max_y < rect2.min_y) | \ - (rect1.min_y > rect2.max_y) - - @staticmethod - def find_intersection(rect1: "Rect", rect2: "Rect") -> Union[None, "Rect"]: - """ - Find the intersection area of two intensity tables and return as Rect object. - If no overlap return none. - """ - if Rect.no_overlap(rect1, rect2): - return None - return Rect(max(rect1.min_x, rect2.min_x), - min(rect1.max_x, rect2.max_x), - max(rect1.min_y, rect2.min_y), - min(rect1.max_y, rect2.max_y)) - - -def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] - ) -> List[Tuple[int, int, "Rect"]]: - """ - Takes a list of xarrays and returns a list of xarray pairs that overlap physically - and their areas of overlap defined by a Rectangle. - """ - all_overlaps: List[Tuple[int, int, "Rect"]] = list() - for idx1, idx2 in itertools.combinations(range(len(xarrays)), 2): - xr1 = xarrays[idx1] - xr2 = xarrays[idx2] - rect1 = Rect(min(xr1['xc']).data, - max(xr1['xc']).data, - min(xr1['yc']).data, - max(xr1['yc']).data) - - rect2 = Rect(min(xr2['xc']).data, - max(xr2['xc']).data, - min(xr2['yc']).data, - max(xr2['yc']).data) - intersection = Rect.find_intersection(rect1, rect2) - if intersection: - all_overlaps.append((idx1, idx2, intersection)) - return all_overlaps - - -def remove_area_of_xarray(it: xr.DataArray, area: Rect) -> xr.DataArray: - """Remove section of xarray within a defined area""" - return it.where((it.xc <= area.min_x) | - (it.xc >= area.max_x) | - (it.yc <= area.min_y) | - (it.yc >= area.max_y), drop=True) - - -def sel_area_of_xarray(it: xr.DataArray, area: Rect) -> xr.DataArray: - """Select on xarray within a defined area""" - return it.where((it.xc > area.min_x) & - (it.xc < area.max_x) & - (it.yc > area.min_y) & - (it.yc < area.max_y), drop=True) - - -def take_max(intersection_rect: Rect, - it1: xr.DataArray, - intersect1: xr.DataArray, - it2: xr.DataArray, - intersect2: xr.DataArray): - # # compare to see which section has more spots - if intersect1.sizes[Features.AXIS] >= intersect2.sizes[Features.AXIS]: - # I1 has more spots remove intesection2 from it2 - it2 = remove_area_of_xarray(it2, intersection_rect) - else: - it1 = remove_area_of_xarray(it1, intersection_rect) - return it1, it2 - - -OVERLAP_STRATEGY_MAP = { - OverlapStrategy.TAKE_MAX: take_max -} - - - - diff --git a/starfish/util/overlap_utils.py b/starfish/util/overlap_utils.py new file mode 100644 index 000000000..c9482ac7a --- /dev/null +++ b/starfish/util/overlap_utils.py @@ -0,0 +1,104 @@ +import itertools +from typing import List, Union, Tuple +import xarray as xr + +from starfish.types import Coordinates, Features +from starfish.types._constants import OverlapStrategy + + +class Area: + def __init__(self, min_x, max_x, min_y, max_y): + """ + Small class that defines an area of physical space by its bottom left + and top right coordinates. + """ + self.min_x = min_x + self.max_x = max_x + self.min_y = min_y + self.max_y = max_y + + @staticmethod + def no_overlap(area1: "Area", area2: "Area") -> bool: + """Return True if two rectangles do not overlap""" + return (area1.max_x < area2.min_x) | \ + (area1.min_x > area2.max_x) | \ + (area1.max_y < area2.min_y) | \ + (area1.min_y > area2.max_y) + + @staticmethod + def find_intersection(area1: "Area", area2: "Area") -> Union[None, "Area"]: + """ + Find the intersection area of two areas and return as new Area object. + If no overlap return none. + """ + if Area.no_overlap(area1, area2): + return None + return Area(max(area1.min_x, area2.min_x), + min(area1.max_x, area2.max_x), + max(area1.min_y, area2.min_y), + min(area1.max_y, area2.max_y)) + + +def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] + ) -> List[Tuple[int, int, "Area"]]: + """ + Takes a list of xarrays and returns a list of xarray index pairs that overlap physically + and their area of overlap defined by an Area object. + """ + all_overlaps: List[Tuple[int, int, "Area"]] = list() + for idx1, idx2 in itertools.combinations(range(len(xarrays)), 2): + xr1 = xarrays[idx1] + xr2 = xarrays[idx2] + area1 = Area(min(xr1[Coordinates.X.value]).data, + max(xr1[Coordinates.X.value]).data, + min(xr1[Coordinates.Y.value]).data, + max(xr1[Coordinates.Y.value]).data) + + area2 = Area(min(xr2[Coordinates.X.value]).data, + max(xr2[Coordinates.X.value]).data, + min(xr2[Coordinates.Y.value]).data, + max(xr2[Coordinates.Y.value]).data) + intersection = Area.find_intersection(area1, area2) + if intersection: + all_overlaps.append((idx1, idx2, intersection)) + return all_overlaps + + +def remove_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: + """Remove area from xarray""" + return it.where((it.xc <= area.min_x) | + (it.xc >= area.max_x) | + (it.yc <= area.min_y) | + (it.yc >= area.max_y), drop=True) + + +def sel_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: + """Select on xarray within a defined area""" + return it.where((it.xc > area.min_x) & + (it.xc < area.max_x) & + (it.yc > area.min_y) & + (it.yc < area.max_y), drop=True) + + +def take_max(intersection_rect: Area, + it1: xr.DataArray, + it2: xr.DataArray + ): + # # compare to see which section has more spots + intersect1 = sel_area_of_xarray(it1, intersection_rect) + intersect2 = sel_area_of_xarray(it2, intersection_rect) + if intersect1.sizes[Features.AXIS] >= intersect2.sizes[Features.AXIS]: + # I1 has more spots remove intesection2 from it2 + it2 = remove_area_of_xarray(it2, intersection_rect) + else: + it1 = remove_area_of_xarray(it1, intersection_rect) + return it1, it2 + + +OVERLAP_STRATEGY_MAP = { + OverlapStrategy.TAKE_MAX: take_max +} + + + + From c348597f8eeb8fdc957cccaaa9a9f0f3616e9228 Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Wed, 27 Mar 2019 17:01:30 -0700 Subject: [PATCH 15/21] adding better docs, need better tests --- ..._table_concat.py => test_overlap_utils.py} | 21 +++++- starfish/types/_constants.py | 7 +- starfish/util/overlap_utils.py | 71 ++++++++++++++++--- 3 files changed, 85 insertions(+), 14 deletions(-) rename starfish/test/{intensity_table/test_intensity_table_concat.py => test_overlap_utils.py} (83%) diff --git a/starfish/test/intensity_table/test_intensity_table_concat.py b/starfish/test/test_overlap_utils.py similarity index 83% rename from starfish/test/intensity_table/test_intensity_table_concat.py rename to starfish/test/test_overlap_utils.py index 73ca8a58d..49d602daa 100644 --- a/starfish/test/intensity_table/test_intensity_table_concat.py +++ b/starfish/test/test_overlap_utils.py @@ -7,7 +7,23 @@ from starfish.types._constants import OverlapStrategy -def test_overlap(): +def test_find_intersection(): + return None + + +def test_test_no_overlap(): + return None + + +def test_remove_area_of_xarray(): + return None + + +def test_sel_area_of_xarray(): + return None + + +def test_take_max(): codebook = test_utils.codebook_array_factory() it1 = IntensityTable.synthetic_intensities( codebook, @@ -35,5 +51,4 @@ def test_overlap(): process_overlaps=True, overlap_strategy=OverlapStrategy.TAKE_MAX) - assert concatenated.sizes[Features.AXIS] == 26 - + assert concatenated.sizes[Features.AXIS] == 26 \ No newline at end of file diff --git a/starfish/types/_constants.py b/starfish/types/_constants.py index 4315d2f8b..e106bef16 100644 --- a/starfish/types/_constants.py +++ b/starfish/types/_constants.py @@ -78,11 +78,14 @@ class Features: INTENSITY = 'intensity' -class OverlapStrategy: +class OverlapStrategy(AugmentedEnum): + """ + contains options to use when processes physically overlapping IntensityTables + or ImageStacks + """ TAKE_MAX = 'take_max' - class Clip(AugmentedEnum): """ contains clipping options that determine how out-of-bounds values produced by filters are diff --git a/starfish/util/overlap_utils.py b/starfish/util/overlap_utils.py index c9482ac7a..7a86a88e1 100644 --- a/starfish/util/overlap_utils.py +++ b/starfish/util/overlap_utils.py @@ -7,11 +7,11 @@ class Area: + """ + Small class that defines rectangular area of physical space by + its bottom left and top right coordinates. + """ def __init__(self, min_x, max_x, min_y, max_y): - """ - Small class that defines an area of physical space by its bottom left - and top right coordinates. - """ self.min_x = min_x self.max_x = max_x self.min_y = min_y @@ -28,7 +28,7 @@ def no_overlap(area1: "Area", area2: "Area") -> bool: @staticmethod def find_intersection(area1: "Area", area2: "Area") -> Union[None, "Area"]: """ - Find the intersection area of two areas and return as new Area object. + Find the overlap area of two rectangles and return as new Area object. If no overlap return none. """ if Area.no_overlap(area1, area2): @@ -42,8 +42,19 @@ def find_intersection(area1: "Area", area2: "Area") -> Union[None, "Area"]: def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] ) -> List[Tuple[int, int, "Area"]]: """ - Takes a list of xarrays and returns a list of xarray index pairs that overlap physically - and their area of overlap defined by an Area object. + Find all the overlap areas within a list of xarrays. + + Parameters + ---------- + xarrays : List[xr.DataArray] + The list of xarrays to find overlaps in. + + Returns + ------- + List[Tuple[int, int, "Area"]] : + A list of tuples containing the indices of two overlapping + IntensityTables and their Area of intersection. + """ all_overlaps: List[Tuple[int, int, "Area"]] = list() for idx1, idx2 in itertools.combinations(range(len(xarrays)), 2): @@ -65,7 +76,20 @@ def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] def remove_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: - """Remove area from xarray""" + """Return everything in the xarray defined OUTSIDE the input area + + Parameters + ---------- + it: xr.DataArray + The xarray to modify + area: Area + The area to not include in the modified xarray + + Returns + ------- + xr.DataArray : + The xarray without the defined area. + """ return it.where((it.xc <= area.min_x) | (it.xc >= area.max_x) | (it.yc <= area.min_y) | @@ -73,7 +97,20 @@ def remove_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: def sel_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: - """Select on xarray within a defined area""" + """Return everything in the xarray defined WITHIN the input area + + Parameters + ---------- + it: xr.DataArray + The xarray to modify + area: Area + The area to include in the modified xarray + + Returns + ------- + xr.DataArray : + The xarray within the defined area. + """ return it.where((it.xc > area.min_x) & (it.xc < area.max_x) & (it.yc > area.min_y) & @@ -84,6 +121,19 @@ def take_max(intersection_rect: Area, it1: xr.DataArray, it2: xr.DataArray ): + """ + Compare two overlapping xarrays and remove spots from whichever + has less int the overlapping section. + + Parameters + ---------- + intersection_rect : Area + The area of physical overalp between two xarrays + it1 : xr.DataArray + The firs overlapping xarray + it2 : xr.DataArray + The second overlapping xarray + """ # # compare to see which section has more spots intersect1 = sel_area_of_xarray(it1, intersection_rect) intersect2 = sel_area_of_xarray(it2, intersection_rect) @@ -95,6 +145,9 @@ def take_max(intersection_rect: Area, return it1, it2 +""" +The mapping between OverlapStrategy type and the method to use for each. +""" OVERLAP_STRATEGY_MAP = { OverlapStrategy.TAKE_MAX: take_max } From be54805c55042318edc821c4abb831196fc8a781 Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Fri, 29 Mar 2019 15:04:54 -0700 Subject: [PATCH 16/21] added tests need to make work for negative coords --- starfish/intensity_table/intensity_table.py | 3 +- starfish/test/test_overlap_utils.py | 112 +++++++++++++++----- starfish/util/overlap_utils.py | 27 +++-- 3 files changed, 102 insertions(+), 40 deletions(-) diff --git a/starfish/intensity_table/intensity_table.py b/starfish/intensity_table/intensity_table.py index a5ba23de6..4906ec4df 100644 --- a/starfish/intensity_table/intensity_table.py +++ b/starfish/intensity_table/intensity_table.py @@ -402,9 +402,10 @@ def process_overlaps(its: List["IntensityTable"], and process them according to the give overlap strategy """ overlap_pairs = find_overlaps_of_xarrays(its) - for idx1, idx2, intersection_rect in overlap_pairs: + for indexes, intersection_rect in overlap_pairs.items(): overlap_method = OVERLAP_STRATEGY_MAP[overlap_strategy] # modify intensity tables based on overlap strategy + idx1, idx2 = indexes it1, it2 = overlap_method(intersection_rect, its[idx1], its[idx2]) its[idx1] = it1 its[idx2] = it2 diff --git a/starfish/test/test_overlap_utils.py b/starfish/test/test_overlap_utils.py index 49d602daa..e4954e6c1 100644 --- a/starfish/test/test_overlap_utils.py +++ b/starfish/test/test_overlap_utils.py @@ -5,50 +5,104 @@ from starfish.test import test_utils from starfish.types import Coordinates, Features from starfish.types._constants import OverlapStrategy +from starfish.util.overlap_utils import Area, remove_area_of_xarray, sel_area_of_xarray, find_overlaps_of_xarrays -def test_find_intersection(): - return None +def create_intenisty_table_with_coords(area: Area, n_spots=10): + """ + Creates a 50X50 intensity table with physical coordinates within + the given Area. + """ + codebook = test_utils.codebook_array_factory() + it = IntensityTable.synthetic_intensities( + codebook, + num_z=1, + height=50, + width=50, + n_spots=n_spots + ) + # intensity table 1 has 10 spots, xmin = 0, ymin = 0, xmax = 2, ymax = 1 + it[Coordinates.X.value] = xr.DataArray(np.linspace(area.min_x, area.max_x, n_spots), dims=Features.AXIS) + it[Coordinates.Y.value] = xr.DataArray(np.linspace(area.min_y, area.max_y, n_spots), dims=Features.AXIS) + return it -def test_test_no_overlap(): - return None +def test_find_area_intersection(): + area1 = Area(min_x=0, max_x=2, min_y=0, max_y=2) + area2 = Area(min_x=1, max_x=2, min_y=1, max_y=3) + intersection = Area.find_intersection(area1, area2) + # intersection should be area with bottom point (1,1) and top point (2,2) + assert intersection == Area(min_x=1, max_x=2, min_y=1, max_y=2) + + area2 = Area(min_x=3, max_x=5, min_y=3, max_y=5) + intersection = Area.find_intersection(area1, area2) + # no intersection + assert intersection is None + + area2 = Area(min_x=0, max_x=5, min_y=3, max_y=5) + # area 2 right above area one + assert intersection is None + + +def test_find_overlaps_of_xarrays(): + # Create some overlapping intensity tables + it0 = create_intenisty_table_with_coords(Area(min_x=0, max_x=1, min_y=0, max_y=1)) + it1 = create_intenisty_table_with_coords(Area(min_x=.5, max_x=2, min_y=.5, max_y=1.5)) + it2 = create_intenisty_table_with_coords(Area(min_x=1.5, max_x=2.5, min_y=0, max_y=1)) + it3 = create_intenisty_table_with_coords(Area(min_x=0, max_x=1, min_y=1, max_y=2)) + overlaps= find_overlaps_of_xarrays([it0, it1, it2, it3]) + # should have 4 total overlaps + assert len(overlaps) == 4 + # overlap 1 between it0 and it1: + overlap_1 = Area(min_x=.5, max_x=1, min_y=.5, max_y=1) + assert overlap_1 == overlaps[(0, 1)] + # overlap 2 between it1 and it2 + overlap_2 = Area(min_x=1.5, max_x=2, min_y=.5, max_y=1) + assert overlap_2 == overlaps[(1, 2)] + # overlap 3 between it1 and it3 + overlap_3 = Area(min_x=.5, max_x=1, min_y=1, max_y=1.5) + assert overlap_3 == overlaps[(1, 3)] + # overlap 4 between it0 and it3 + overlap_4 = Area(min_x=0, max_x=1, min_y=1, max_y=1) + assert overlap_4 == overlaps[(0, 3)] def test_remove_area_of_xarray(): - return None + it = create_intenisty_table_with_coords(Area(min_x=0, max_x=2, min_y=0, max_y=2), n_spots=10) + + area = Area(min_x=1, max_x=2, min_y=1, max_y=3) + # grab some random coord values in this range + removed_x = it.where(it.xc > 1, drop=True)[Coordinates.X.value].data[0] + removed_y = it.where(it.yc > 1, drop=True)[Coordinates.X.value].data[3] + + it = remove_area_of_xarray(it, area) + # assert coords from removed section are no longer in it + assert not np.any(np.isclose(it[Coordinates.X.value], removed_x)) + assert not np.any(np.isclose(it[Coordinates.Y.value], removed_y)) def test_sel_area_of_xarray(): - return None + it = create_intenisty_table_with_coords(Area(min_x=0, max_x=2, min_y=0, max_y=2), n_spots=10) + area = Area(min_x=1, max_x=2, min_y=1, max_y=3) + it = sel_area_of_xarray(it, area) -def test_take_max(): - codebook = test_utils.codebook_array_factory() - it1 = IntensityTable.synthetic_intensities( - codebook, - num_z=1, - height=50, - width=50, - n_spots=10 - ) - # intensity table 1 has 10 spots, xmin = 0, ymin = 0, xmax = 2, ymax = 1 - it1[Coordinates.X.value] = xr.DataArray(np.linspace(0, 2, 10), dims='features') - it1[Coordinates.Y.value] = xr.DataArray(np.linspace(0, 2, 10), dims='features') + # Assert noew min/max values + assert min(it[Coordinates.X.value]).data >= 1 + assert max(it[Coordinates.X.value]).data <= 2 + assert min(it[Coordinates.Y.value]).data >= 1 + assert max(it[Coordinates.X.value]).data <= 2 - it2 = IntensityTable.synthetic_intensities( - codebook, - num_z=1, - height=50, - width=50, - n_spots=20 - ) - # intensity table 2 has 20 spots, xmin = 1, ymin = 1, xmax = 3, ymax = 3 - it2[Coordinates.X.value] = xr.DataArray(np.linspace(1, 3, 20), dims='features') - it2[Coordinates.Y.value] = xr.DataArray(np.linspace(1, 3, 20), dims='features') + +def test_take_max(): + it1 = create_intenisty_table_with_coords( Area(min_x=0, max_x=2, min_y=0, max_y=2), n_spots=10) + it2 = create_intenisty_table_with_coords(Area(min_x=1, max_x=2, min_y=1, max_y=3), n_spots=20) concatenated = IntensityTable.concatanate_intensity_tables([it1, it2], process_overlaps=True, overlap_strategy=OverlapStrategy.TAKE_MAX) - assert concatenated.sizes[Features.AXIS] == 26 \ No newline at end of file + # The overlap section hits half of the spots from each intensity table, 5 from it1 and 10 from i21. + # It2 wins and the resulting concatenated table should have all the spots from it2 (20) and 6 from + # it1 (6) for a total of 26 spots + assert concatenated.sizes[Features.AXIS] == 26 diff --git a/starfish/util/overlap_utils.py b/starfish/util/overlap_utils.py index 7a86a88e1..d7d744a43 100644 --- a/starfish/util/overlap_utils.py +++ b/starfish/util/overlap_utils.py @@ -1,5 +1,5 @@ import itertools -from typing import List, Union, Tuple +from typing import List, Dict, Union, Tuple import xarray as xr from starfish.types import Coordinates, Features @@ -17,9 +17,16 @@ def __init__(self, min_x, max_x, min_y, max_y): self.min_y = min_y self.max_y = max_y + def __eq__(self, other): + return self.min_x == other.min_x and \ + self.min_y == other.min_y and \ + self.max_x == other.max_x and\ + self.max_y == other.max_y + @staticmethod - def no_overlap(area1: "Area", area2: "Area") -> bool: + def _no_overlap(area1: "Area", area2: "Area") -> bool: """Return True if two rectangles do not overlap""" + # TODO ACCOUNT FOR NEGATIVES return (area1.max_x < area2.min_x) | \ (area1.min_x > area2.max_x) | \ (area1.max_y < area2.min_y) | \ @@ -31,7 +38,7 @@ def find_intersection(area1: "Area", area2: "Area") -> Union[None, "Area"]: Find the overlap area of two rectangles and return as new Area object. If no overlap return none. """ - if Area.no_overlap(area1, area2): + if Area._no_overlap(area1, area2): return None return Area(max(area1.min_x, area2.min_x), min(area1.max_x, area2.max_x), @@ -40,7 +47,7 @@ def find_intersection(area1: "Area", area2: "Area") -> Union[None, "Area"]: def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] - ) -> List[Tuple[int, int, "Area"]]: + ) -> Dict[Tuple[int, int], "Area"]: """ Find all the overlap areas within a list of xarrays. @@ -56,7 +63,7 @@ def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] IntensityTables and their Area of intersection. """ - all_overlaps: List[Tuple[int, int, "Area"]] = list() + all_overlaps: Dict[Tuple[int, int], "Area"] = dict() for idx1, idx2 in itertools.combinations(range(len(xarrays)), 2): xr1 = xarrays[idx1] xr2 = xarrays[idx2] @@ -71,7 +78,7 @@ def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] max(xr2[Coordinates.Y.value]).data) intersection = Area.find_intersection(area1, area2) if intersection: - all_overlaps.append((idx1, idx2, intersection)) + all_overlaps[(idx1, idx2)] = intersection return all_overlaps @@ -111,10 +118,10 @@ def sel_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: xr.DataArray : The xarray within the defined area. """ - return it.where((it.xc > area.min_x) & - (it.xc < area.max_x) & - (it.yc > area.min_y) & - (it.yc < area.max_y), drop=True) + return it.where((it.xc >= area.min_x) & + (it.xc <= area.max_x) & + (it.yc >= area.min_y) & + (it.yc <= area.max_y), drop=True) def take_max(intersection_rect: Area, From 97fdc3554380aa3144bf611399c116d07fd15f45 Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Fri, 29 Mar 2019 16:25:22 -0700 Subject: [PATCH 17/21] fixing lint errors --- starfish/intensity_table/intensity_table.py | 8 +- starfish/test/test_overlap_utils.py | 83 ++++++++++++++++----- starfish/types/_constants.py | 1 - starfish/util/overlap_utils.py | 57 +++++++------- 4 files changed, 97 insertions(+), 52 deletions(-) diff --git a/starfish/intensity_table/intensity_table.py b/starfish/intensity_table/intensity_table.py index 4906ec4df..ca3e682c4 100644 --- a/starfish/intensity_table/intensity_table.py +++ b/starfish/intensity_table/intensity_table.py @@ -8,11 +8,10 @@ import xarray as xr from starfish.expression_matrix.expression_matrix import ExpressionMatrix - -from starfish.types._constants import OverlapStrategy from starfish.types import Axes, DecodedSpots, Features, LOG, SpotAttributes, STARFISH_EXTRAS_KEY +from starfish.types._constants import OverlapStrategy from starfish.util.dtype import preserve_float_range -from starfish.util.overlap_utils import( +from starfish.util.overlap_utils import ( find_overlaps_of_xarrays, OVERLAP_STRATEGY_MAP, ) @@ -413,9 +412,8 @@ def process_overlaps(its: List["IntensityTable"], @staticmethod def concatanate_intensity_tables(intensity_tables: List["IntensityTable"], - process_overlaps: bool = False, overlap_strategy: OverlapStrategy = None): - if process_overlaps: + if overlap_strategy: intensity_tables = IntensityTable.process_overlaps(intensity_tables, overlap_strategy) return xr.concat(intensity_tables, dim=Features.AXIS) diff --git a/starfish/test/test_overlap_utils.py b/starfish/test/test_overlap_utils.py index e4954e6c1..720709fca 100644 --- a/starfish/test/test_overlap_utils.py +++ b/starfish/test/test_overlap_utils.py @@ -5,13 +5,25 @@ from starfish.test import test_utils from starfish.types import Coordinates, Features from starfish.types._constants import OverlapStrategy -from starfish.util.overlap_utils import Area, remove_area_of_xarray, sel_area_of_xarray, find_overlaps_of_xarrays +from starfish.util.overlap_utils import ( + Area, + find_overlaps_of_xarrays, + remove_area_of_xarray, + sel_area_of_xarray +) def create_intenisty_table_with_coords(area: Area, n_spots=10): """ Creates a 50X50 intensity table with physical coordinates within the given Area. + + Parameters + ---------- + area: Area + The area of physical space the IntensityTable should be defined over + n_spots: + Number of spots to add to the IntensityTable """ codebook = test_utils.codebook_array_factory() it = IntensityTable.synthetic_intensities( @@ -22,12 +34,17 @@ def create_intenisty_table_with_coords(area: Area, n_spots=10): n_spots=n_spots ) # intensity table 1 has 10 spots, xmin = 0, ymin = 0, xmax = 2, ymax = 1 - it[Coordinates.X.value] = xr.DataArray(np.linspace(area.min_x, area.max_x, n_spots), dims=Features.AXIS) - it[Coordinates.Y.value] = xr.DataArray(np.linspace(area.min_y, area.max_y, n_spots), dims=Features.AXIS) + it[Coordinates.X.value] = xr.DataArray(np.linspace(area.min_x, area.max_x, n_spots), + dims=Features.AXIS) + it[Coordinates.Y.value] = xr.DataArray(np.linspace(area.min_y, area.max_y, n_spots), + dims=Features.AXIS) return it def test_find_area_intersection(): + """ + Create various Area objects and verify their intersection are calculated correctly + """ area1 = Area(min_x=0, max_x=2, min_y=0, max_y=2) area2 = Area(min_x=1, max_x=2, min_y=1, max_y=3) intersection = Area.find_intersection(area1, area2) @@ -40,17 +57,36 @@ def test_find_area_intersection(): assert intersection is None area2 = Area(min_x=0, max_x=5, min_y=3, max_y=5) + intersection = Area.find_intersection(area1, area2) # area 2 right above area one assert intersection is None + # try negatives + area1 = Area(min_x=-1, max_x=1, min_y=0, max_y=2) + area2 = Area(min_x=0, max_x=2, min_y=0, max_y=2) + intersection = Area.find_intersection(area1, area2) + assert intersection == Area(min_x=0, max_x=1, min_y=0, max_y=2) + + area2 = Area(min_x=-3, max_x=-2, min_y=0, max_y=2) + intersection = Area.find_intersection(area1, area2) + assert intersection is None + def test_find_overlaps_of_xarrays(): + """ + Create a list of overlapping IntensityTables and verify we identify the correct + overlapping sections + """ # Create some overlapping intensity tables - it0 = create_intenisty_table_with_coords(Area(min_x=0, max_x=1, min_y=0, max_y=1)) - it1 = create_intenisty_table_with_coords(Area(min_x=.5, max_x=2, min_y=.5, max_y=1.5)) - it2 = create_intenisty_table_with_coords(Area(min_x=1.5, max_x=2.5, min_y=0, max_y=1)) - it3 = create_intenisty_table_with_coords(Area(min_x=0, max_x=1, min_y=1, max_y=2)) - overlaps= find_overlaps_of_xarrays([it0, it1, it2, it3]) + it0 = create_intenisty_table_with_coords(Area(min_x=0, max_x=1, + min_y=0, max_y=1)) + it1 = create_intenisty_table_with_coords(Area(min_x=.5, max_x=2, + min_y=.5, max_y=1.5)) + it2 = create_intenisty_table_with_coords(Area(min_x=1.5, max_x=2.5, + min_y=0, max_y=1)) + it3 = create_intenisty_table_with_coords(Area(min_x=0, max_x=1, + min_y=1, max_y=2)) + overlaps = find_overlaps_of_xarrays([it0, it1, it2, it3]) # should have 4 total overlaps assert len(overlaps) == 4 # overlap 1 between it0 and it1: @@ -68,7 +104,11 @@ def test_find_overlaps_of_xarrays(): def test_remove_area_of_xarray(): - it = create_intenisty_table_with_coords(Area(min_x=0, max_x=2, min_y=0, max_y=2), n_spots=10) + """ + Tests removing a section of an IntensityTable defined by its physical area + """ + it = create_intenisty_table_with_coords(Area(min_x=0, max_x=2, + min_y=0, max_y=2), n_spots=10) area = Area(min_x=1, max_x=2, min_y=1, max_y=3) # grab some random coord values in this range @@ -82,6 +122,9 @@ def test_remove_area_of_xarray(): def test_sel_area_of_xarray(): + """ + Tests selecting a section of an IntensityTable defined by its physical area + """ it = create_intenisty_table_with_coords(Area(min_x=0, max_x=2, min_y=0, max_y=2), n_spots=10) area = Area(min_x=1, max_x=2, min_y=1, max_y=3) @@ -95,14 +138,20 @@ def test_sel_area_of_xarray(): def test_take_max(): - it1 = create_intenisty_table_with_coords( Area(min_x=0, max_x=2, min_y=0, max_y=2), n_spots=10) - it2 = create_intenisty_table_with_coords(Area(min_x=1, max_x=2, min_y=1, max_y=3), n_spots=20) + """ + Create two overlapping IntensityTables with differing number of spots and verify that + by concatenating them with the TAKE_MAX strategy we only include spots in the overlapping + section from the IntensityTable that had the most. + """ + it1 = create_intenisty_table_with_coords(Area(min_x=0, max_x=2, + min_y=0, max_y=2), n_spots=10) + it2 = create_intenisty_table_with_coords(Area(min_x=1, max_x=2, + min_y=1, max_y=3), n_spots=20) - concatenated = IntensityTable.concatanate_intensity_tables([it1, it2], - process_overlaps=True, - overlap_strategy=OverlapStrategy.TAKE_MAX) + concatenated = IntensityTable.concatanate_intensity_tables( + [it1, it2], overlap_strategy=OverlapStrategy.TAKE_MAX) - # The overlap section hits half of the spots from each intensity table, 5 from it1 and 10 from i21. - # It2 wins and the resulting concatenated table should have all the spots from it2 (20) and 6 from - # it1 (6) for a total of 26 spots + # The overlap section hits half of the spots from each intensity table, 5 from it1 + # and 10 from i21. It2 wins and the resulting concatenated table should have all the + # spots from it2 (20) and 6 from it1 (6) for a total of 26 spots assert concatenated.sizes[Features.AXIS] == 26 diff --git a/starfish/types/_constants.py b/starfish/types/_constants.py index e106bef16..c189aef57 100644 --- a/starfish/types/_constants.py +++ b/starfish/types/_constants.py @@ -94,4 +94,3 @@ class Clip(AugmentedEnum): CLIP = 'clip' SCALE_BY_IMAGE = 'scale_by_image' SCALE_BY_CHUNK = 'scale_by_chunk' - diff --git a/starfish/util/overlap_utils.py b/starfish/util/overlap_utils.py index d7d744a43..b4180db02 100644 --- a/starfish/util/overlap_utils.py +++ b/starfish/util/overlap_utils.py @@ -1,5 +1,6 @@ import itertools -from typing import List, Dict, Union, Tuple +from typing import Dict, List, Tuple, Union + import xarray as xr from starfish.types import Coordinates, Features @@ -18,19 +19,20 @@ def __init__(self, min_x, max_x, min_y, max_y): self.max_y = max_y def __eq__(self, other): - return self.min_x == other.min_x and \ - self.min_y == other.min_y and \ - self.max_x == other.max_x and\ - self.max_y == other.max_y + return (self.min_x == other.min_x + and self.min_y == other.min_y + and self.max_x == other.max_x + and self.max_y == other.max_y) @staticmethod - def _no_overlap(area1: "Area", area2: "Area") -> bool: + def _overlap(area1: "Area", area2: "Area") -> bool: """Return True if two rectangles do not overlap""" # TODO ACCOUNT FOR NEGATIVES - return (area1.max_x < area2.min_x) | \ - (area1.min_x > area2.max_x) | \ - (area1.max_y < area2.min_y) | \ - (area1.min_y > area2.max_y) + if (area1.max_x < area2.min_x) or (area1.min_x > area2.max_x): + return False + if (area1.max_y < area2.min_y) or (area1.min_y > area2.max_y): + return False + return True @staticmethod def find_intersection(area1: "Area", area2: "Area") -> Union[None, "Area"]: @@ -38,12 +40,12 @@ def find_intersection(area1: "Area", area2: "Area") -> Union[None, "Area"]: Find the overlap area of two rectangles and return as new Area object. If no overlap return none. """ - if Area._no_overlap(area1, area2): - return None - return Area(max(area1.min_x, area2.min_x), - min(area1.max_x, area2.max_x), - max(area1.min_y, area2.min_y), - min(area1.max_y, area2.max_y)) + if Area._overlap(area1, area2): + return Area(max(area1.min_x, area2.min_x), + min(area1.max_x, area2.max_x), + max(area1.min_y, area2.min_y), + min(area1.max_y, area2.max_y)) + return None def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] @@ -63,7 +65,7 @@ def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] IntensityTables and their Area of intersection. """ - all_overlaps: Dict[Tuple[int, int], "Area"] = dict() + all_overlaps: Dict[Tuple[int, int], "Area"] = dict() for idx1, idx2 in itertools.combinations(range(len(xarrays)), 2): xr1 = xarrays[idx1] xr2 = xarrays[idx2] @@ -97,10 +99,11 @@ def remove_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: xr.DataArray : The xarray without the defined area. """ - return it.where((it.xc <= area.min_x) | - (it.xc >= area.max_x) | - (it.yc <= area.min_y) | - (it.yc >= area.max_y), drop=True) + return it.where((it.xc <= area.min_x) + | (it.xc >= area.max_x) + | (it.yc <= area.min_y) + | (it.yc >= area.max_y), + drop=True) def sel_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: @@ -118,10 +121,10 @@ def sel_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: xr.DataArray : The xarray within the defined area. """ - return it.where((it.xc >= area.min_x) & - (it.xc <= area.max_x) & - (it.yc >= area.min_y) & - (it.yc <= area.max_y), drop=True) + return it.where((it.xc >= area.min_x) + & (it.xc <= area.max_x) + & (it.yc >= area.min_y) + & (it.yc <= area.max_y), drop=True) def take_max(intersection_rect: Area, @@ -158,7 +161,3 @@ def take_max(intersection_rect: Area, OVERLAP_STRATEGY_MAP = { OverlapStrategy.TAKE_MAX: take_max } - - - - From 6499c7a2d281bed4fc16a6b0f6cf64dbd6e3ea0d Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Fri, 29 Mar 2019 16:27:34 -0700 Subject: [PATCH 18/21] notebooks --- .../DARTFISH_Pipeline_-_Human_Occipital_Cortex_-_1_FOV.ipynb | 2 +- notebooks/MERFISH_Pipeline_-_U2O2_Cell_Culture_-_1_FOV.ipynb | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/notebooks/DARTFISH_Pipeline_-_Human_Occipital_Cortex_-_1_FOV.ipynb b/notebooks/DARTFISH_Pipeline_-_Human_Occipital_Cortex_-_1_FOV.ipynb index 7a99c2132..f0ed993ae 100644 --- a/notebooks/DARTFISH_Pipeline_-_Human_Occipital_Cortex_-_1_FOV.ipynb +++ b/notebooks/DARTFISH_Pipeline_-_Human_Occipital_Cortex_-_1_FOV.ipynb @@ -421,4 +421,4 @@ }, "nbformat": 4, "nbformat_minor": 2 -} \ No newline at end of file +} diff --git a/notebooks/MERFISH_Pipeline_-_U2O2_Cell_Culture_-_1_FOV.ipynb b/notebooks/MERFISH_Pipeline_-_U2O2_Cell_Culture_-_1_FOV.ipynb index ed072eb0c..ec746c6e8 100644 --- a/notebooks/MERFISH_Pipeline_-_U2O2_Cell_Culture_-_1_FOV.ipynb +++ b/notebooks/MERFISH_Pipeline_-_U2O2_Cell_Culture_-_1_FOV.ipynb @@ -402,4 +402,4 @@ }, "nbformat": 4, "nbformat_minor": 2 -} \ No newline at end of file +} From 847ef42c0d032154aadb7fd86efd62727bc8b4fd Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Fri, 29 Mar 2019 16:44:22 -0700 Subject: [PATCH 19/21] small fixes --- starfish/intensity_table/intensity_table.py | 9 +++++---- starfish/test/test_overlap_utils.py | 2 +- starfish/util/overlap_utils.py | 13 ++++++------- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/starfish/intensity_table/intensity_table.py b/starfish/intensity_table/intensity_table.py index ca3e682c4..498ea2ea3 100644 --- a/starfish/intensity_table/intensity_table.py +++ b/starfish/intensity_table/intensity_table.py @@ -401,11 +401,12 @@ def process_overlaps(its: List["IntensityTable"], and process them according to the give overlap strategy """ overlap_pairs = find_overlaps_of_xarrays(its) - for indexes, intersection_rect in overlap_pairs.items(): + for indices, intersection_area in overlap_pairs.items(): overlap_method = OVERLAP_STRATEGY_MAP[overlap_strategy] - # modify intensity tables based on overlap strategy - idx1, idx2 = indexes - it1, it2 = overlap_method(intersection_rect, its[idx1], its[idx2]) + idx1, idx2 = indices + # modify IntensityTables based on overlap strategy + it1, it2 = overlap_method(intersection_area, its[idx1], its[idx2]) + # replace IntensityTables in list its[idx1] = it1 its[idx2] = it2 return its diff --git a/starfish/test/test_overlap_utils.py b/starfish/test/test_overlap_utils.py index 720709fca..d7b9ceb9d 100644 --- a/starfish/test/test_overlap_utils.py +++ b/starfish/test/test_overlap_utils.py @@ -130,7 +130,7 @@ def test_sel_area_of_xarray(): area = Area(min_x=1, max_x=2, min_y=1, max_y=3) it = sel_area_of_xarray(it, area) - # Assert noew min/max values + # Assert new min/max values assert min(it[Coordinates.X.value]).data >= 1 assert max(it[Coordinates.X.value]).data <= 2 assert min(it[Coordinates.Y.value]).data >= 1 diff --git a/starfish/util/overlap_utils.py b/starfish/util/overlap_utils.py index b4180db02..96b9fe41b 100644 --- a/starfish/util/overlap_utils.py +++ b/starfish/util/overlap_utils.py @@ -26,8 +26,7 @@ def __eq__(self, other): @staticmethod def _overlap(area1: "Area", area2: "Area") -> bool: - """Return True if two rectangles do not overlap""" - # TODO ACCOUNT FOR NEGATIVES + """Return True if two rectangles overlap""" if (area1.max_x < area2.min_x) or (area1.min_x > area2.max_x): return False if (area1.max_y < area2.min_y) or (area1.min_y > area2.max_y): @@ -60,8 +59,8 @@ def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] Returns ------- - List[Tuple[int, int, "Area"]] : - A list of tuples containing the indices of two overlapping + Dict[Tuple[int, int], "Area"]] : + A dictionary of tuples containing the indices of two overlapping IntensityTables and their Area of intersection. """ @@ -133,14 +132,14 @@ def take_max(intersection_rect: Area, ): """ Compare two overlapping xarrays and remove spots from whichever - has less int the overlapping section. + has less in the overlapping section. Parameters ---------- intersection_rect : Area - The area of physical overalp between two xarrays + The area of physical overalap between two xarrays it1 : xr.DataArray - The firs overlapping xarray + The first overlapping xarray it2 : xr.DataArray The second overlapping xarray """ From be5b9e184291cd26d82dc38ca36e8b2f243f1cef Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Wed, 3 Apr 2019 09:21:11 -0700 Subject: [PATCH 20/21] review changes --- starfish/intensity_table/intensity_table.py | 8 +-- starfish/test/test_overlap_utils.py | 34 ++++++------- starfish/util/overlap_utils.py | 55 ++++++++++++--------- 3 files changed, 50 insertions(+), 47 deletions(-) diff --git a/starfish/intensity_table/intensity_table.py b/starfish/intensity_table/intensity_table.py index 498ea2ea3..cebb14a03 100644 --- a/starfish/intensity_table/intensity_table.py +++ b/starfish/intensity_table/intensity_table.py @@ -397,15 +397,15 @@ def from_image_stack( def process_overlaps(its: List["IntensityTable"], overlap_strategy: OverlapStrategy ) -> List["IntensityTable"]: - """Find the overlapping sections between IntensityTables - and process them according to the give overlap strategy + """Find the overlapping sections between IntensityTables and process them according + to the given overlap strategy """ overlap_pairs = find_overlaps_of_xarrays(its) - for indices, intersection_area in overlap_pairs.items(): + for indices in overlap_pairs: overlap_method = OVERLAP_STRATEGY_MAP[overlap_strategy] idx1, idx2 = indices # modify IntensityTables based on overlap strategy - it1, it2 = overlap_method(intersection_area, its[idx1], its[idx2]) + it1, it2 = overlap_method(its[idx1], its[idx2]) # replace IntensityTables in list its[idx1] = it1 its[idx2] = it2 diff --git a/starfish/test/test_overlap_utils.py b/starfish/test/test_overlap_utils.py index d7b9ceb9d..be35a1474 100644 --- a/starfish/test/test_overlap_utils.py +++ b/starfish/test/test_overlap_utils.py @@ -13,7 +13,7 @@ ) -def create_intenisty_table_with_coords(area: Area, n_spots=10): +def create_intensity_table_with_coords(area: Area, n_spots: int=10) -> IntensityTable: """ Creates a 50X50 intensity table with physical coordinates within the given Area. @@ -78,36 +78,32 @@ def test_find_overlaps_of_xarrays(): overlapping sections """ # Create some overlapping intensity tables - it0 = create_intenisty_table_with_coords(Area(min_x=0, max_x=1, + it0 = create_intensity_table_with_coords(Area(min_x=0, max_x=1, min_y=0, max_y=1)) - it1 = create_intenisty_table_with_coords(Area(min_x=.5, max_x=2, + it1 = create_intensity_table_with_coords(Area(min_x=.5, max_x=2, min_y=.5, max_y=1.5)) - it2 = create_intenisty_table_with_coords(Area(min_x=1.5, max_x=2.5, + it2 = create_intensity_table_with_coords(Area(min_x=1.5, max_x=2.5, min_y=0, max_y=1)) - it3 = create_intenisty_table_with_coords(Area(min_x=0, max_x=1, + it3 = create_intensity_table_with_coords(Area(min_x=0, max_x=1, min_y=1, max_y=2)) overlaps = find_overlaps_of_xarrays([it0, it1, it2, it3]) # should have 4 total overlaps assert len(overlaps) == 4 # overlap 1 between it0 and it1: - overlap_1 = Area(min_x=.5, max_x=1, min_y=.5, max_y=1) - assert overlap_1 == overlaps[(0, 1)] - # overlap 2 between it1 and it2 - overlap_2 = Area(min_x=1.5, max_x=2, min_y=.5, max_y=1) - assert overlap_2 == overlaps[(1, 2)] + assert (0, 1) in overlaps + # overlap 1 between it0 and it1: + assert (1, 2) in overlaps # overlap 3 between it1 and it3 - overlap_3 = Area(min_x=.5, max_x=1, min_y=1, max_y=1.5) - assert overlap_3 == overlaps[(1, 3)] + assert (1, 3) in overlaps # overlap 4 between it0 and it3 - overlap_4 = Area(min_x=0, max_x=1, min_y=1, max_y=1) - assert overlap_4 == overlaps[(0, 3)] + assert (0, 3) in overlaps def test_remove_area_of_xarray(): """ Tests removing a section of an IntensityTable defined by its physical area """ - it = create_intenisty_table_with_coords(Area(min_x=0, max_x=2, + it = create_intensity_table_with_coords(Area(min_x=0, max_x=2, min_y=0, max_y=2), n_spots=10) area = Area(min_x=1, max_x=2, min_y=1, max_y=3) @@ -125,7 +121,7 @@ def test_sel_area_of_xarray(): """ Tests selecting a section of an IntensityTable defined by its physical area """ - it = create_intenisty_table_with_coords(Area(min_x=0, max_x=2, min_y=0, max_y=2), n_spots=10) + it = create_intensity_table_with_coords(Area(min_x=0, max_x=2, min_y=0, max_y=2), n_spots=10) area = Area(min_x=1, max_x=2, min_y=1, max_y=3) it = sel_area_of_xarray(it, area) @@ -143,9 +139,9 @@ def test_take_max(): by concatenating them with the TAKE_MAX strategy we only include spots in the overlapping section from the IntensityTable that had the most. """ - it1 = create_intenisty_table_with_coords(Area(min_x=0, max_x=2, + it1 = create_intensity_table_with_coords(Area(min_x=0, max_x=2, min_y=0, max_y=2), n_spots=10) - it2 = create_intenisty_table_with_coords(Area(min_x=1, max_x=2, + it2 = create_intensity_table_with_coords(Area(min_x=1, max_x=2, min_y=1, max_y=3), n_spots=20) concatenated = IntensityTable.concatanate_intensity_tables( @@ -153,5 +149,5 @@ def test_take_max(): # The overlap section hits half of the spots from each intensity table, 5 from it1 # and 10 from i21. It2 wins and the resulting concatenated table should have all the - # spots from it2 (20) and 6 from it1 (6) for a total of 26 spots + # spots from it2 (20) and 6 (one on the boarder) from it1 (6) for a total of 26 spots assert concatenated.sizes[Features.AXIS] == 26 diff --git a/starfish/util/overlap_utils.py b/starfish/util/overlap_utils.py index 96b9fe41b..f4e423ff9 100644 --- a/starfish/util/overlap_utils.py +++ b/starfish/util/overlap_utils.py @@ -1,9 +1,9 @@ import itertools -from typing import Dict, List, Tuple, Union +from typing import List, Optional, Tuple import xarray as xr -from starfish.types import Coordinates, Features +from starfish.types import Coordinates, Features, Number from starfish.types._constants import OverlapStrategy @@ -12,7 +12,7 @@ class Area: Small class that defines rectangular area of physical space by its bottom left and top right coordinates. """ - def __init__(self, min_x, max_x, min_y, max_y): + def __init__(self, min_x: Number, max_x: Number, min_y: Number, max_y: Number): self.min_x = min_x self.max_x = max_x self.min_y = min_y @@ -34,10 +34,10 @@ def _overlap(area1: "Area", area2: "Area") -> bool: return True @staticmethod - def find_intersection(area1: "Area", area2: "Area") -> Union[None, "Area"]: + def find_intersection(area1: "Area", area2: "Area") -> Optional["Area"]: """ Find the overlap area of two rectangles and return as new Area object. - If no overlap return none. + If the two rectangles do not overlap, return None. """ if Area._overlap(area1, area2): return Area(max(area1.min_x, area2.min_x), @@ -47,8 +47,7 @@ def find_intersection(area1: "Area", area2: "Area") -> Union[None, "Area"]: return None -def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] - ) -> Dict[Tuple[int, int], "Area"]: +def find_overlaps_of_xarrays(xarrays: List[xr.DataArray]) -> List[Tuple[int, int]]: """ Find all the overlap areas within a list of xarrays. @@ -59,12 +58,12 @@ def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] Returns ------- - Dict[Tuple[int, int], "Area"]] : - A dictionary of tuples containing the indices of two overlapping - IntensityTables and their Area of intersection. + list[Tuple[int, int]] : + A list of tuples containing the indices of two overlapping + IntensityTables. """ - all_overlaps: Dict[Tuple[int, int], "Area"] = dict() + all_overlaps: List[Tuple[int, int]] = list() for idx1, idx2 in itertools.combinations(range(len(xarrays)), 2): xr1 = xarrays[idx1] xr2 = xarrays[idx2] @@ -77,14 +76,14 @@ def find_overlaps_of_xarrays(xarrays: List[xr.DataArray] max(xr2[Coordinates.X.value]).data, min(xr2[Coordinates.Y.value]).data, max(xr2[Coordinates.Y.value]).data) - intersection = Area.find_intersection(area1, area2) - if intersection: - all_overlaps[(idx1, idx2)] = intersection + if Area._overlap(area1, area2): + all_overlaps.append((idx1, idx2)) return all_overlaps def remove_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: """Return everything in the xarray defined OUTSIDE the input area + including values on the boundary of the area defined. Parameters ---------- @@ -95,7 +94,7 @@ def remove_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: Returns ------- - xr.DataArray : + xr.DataArray : The xarray without the defined area. """ return it.where((it.xc <= area.min_x) @@ -107,6 +106,7 @@ def remove_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: def sel_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: """Return everything in the xarray defined WITHIN the input area + including values on the boundary of the area defined. Parameters ---------- @@ -126,26 +126,33 @@ def sel_area_of_xarray(it: xr.DataArray, area: Area) -> xr.DataArray: & (it.yc <= area.max_y), drop=True) -def take_max(intersection_rect: Area, - it1: xr.DataArray, - it2: xr.DataArray - ): +def take_max(it1: xr.DataArray, it2: xr.DataArray): """ - Compare two overlapping xarrays and remove spots from whichever - has less in the overlapping section. + Compare two overlapping xarrays and remove spots from whichever has less in the + overlapping section. Parameters ---------- - intersection_rect : Area - The area of physical overalap between two xarrays it1 : xr.DataArray The first overlapping xarray it2 : xr.DataArray The second overlapping xarray """ - # # compare to see which section has more spots + area1 = Area(min(it1[Coordinates.X.value]).data, + max(it1[Coordinates.X.value]).data, + min(it1[Coordinates.Y.value]).data, + max(it1[Coordinates.Y.value]).data) + + area2 = Area(min(it2[Coordinates.X.value]).data, + max(it2[Coordinates.X.value]).data, + min(it2[Coordinates.Y.value]).data, + max(it2[Coordinates.Y.value]).data) + intersection_rect = Area.find_intersection(area1, area2) + if not intersection_rect: + raise ValueError("The given xarrays do not overlap") intersect1 = sel_area_of_xarray(it1, intersection_rect) intersect2 = sel_area_of_xarray(it2, intersection_rect) + # compare to see which section has more spots if intersect1.sizes[Features.AXIS] >= intersect2.sizes[Features.AXIS]: # I1 has more spots remove intesection2 from it2 it2 = remove_area_of_xarray(it2, intersection_rect) From a7523732ba3f2d80c319d6027ad754e2d309293a Mon Sep 17 00:00:00 2001 From: Shannon Axelrod Date: Wed, 10 Apr 2019 10:33:29 -0700 Subject: [PATCH 21/21] review comments --- starfish/intensity_table/intensity_table.py | 16 ++++++++-------- starfish/test/test_overlap_utils.py | 2 +- starfish/types/__init__.py | 1 + starfish/util/overlap_utils.py | 10 +++++----- 4 files changed, 15 insertions(+), 14 deletions(-) diff --git a/starfish/intensity_table/intensity_table.py b/starfish/intensity_table/intensity_table.py index cebb14a03..cdcb358d7 100644 --- a/starfish/intensity_table/intensity_table.py +++ b/starfish/intensity_table/intensity_table.py @@ -9,7 +9,7 @@ from starfish.expression_matrix.expression_matrix import ExpressionMatrix from starfish.types import Axes, DecodedSpots, Features, LOG, SpotAttributes, STARFISH_EXTRAS_KEY -from starfish.types._constants import OverlapStrategy +from starfish.types import OverlapStrategy from starfish.util.dtype import preserve_float_range from starfish.util.overlap_utils import ( find_overlaps_of_xarrays, @@ -394,26 +394,26 @@ def from_image_stack( return IntensityTable.from_spot_data(intensity_data, pixel_coordinates) @staticmethod - def process_overlaps(its: List["IntensityTable"], + def process_overlaps(intensity_tables: List["IntensityTable"], overlap_strategy: OverlapStrategy ) -> List["IntensityTable"]: """Find the overlapping sections between IntensityTables and process them according to the given overlap strategy """ - overlap_pairs = find_overlaps_of_xarrays(its) + overlap_pairs = find_overlaps_of_xarrays(intensity_tables) for indices in overlap_pairs: overlap_method = OVERLAP_STRATEGY_MAP[overlap_strategy] idx1, idx2 = indices # modify IntensityTables based on overlap strategy - it1, it2 = overlap_method(its[idx1], its[idx2]) + it1, it2 = overlap_method(intensity_tables[idx1], intensity_tables[idx2]) # replace IntensityTables in list - its[idx1] = it1 - its[idx2] = it2 - return its + intensity_tables[idx1] = it1 + intensity_tables[idx2] = it2 + return intensity_tables @staticmethod def concatanate_intensity_tables(intensity_tables: List["IntensityTable"], - overlap_strategy: OverlapStrategy = None): + overlap_strategy: Optional[OverlapStrategy] = None): if overlap_strategy: intensity_tables = IntensityTable.process_overlaps(intensity_tables, overlap_strategy) diff --git a/starfish/test/test_overlap_utils.py b/starfish/test/test_overlap_utils.py index be35a1474..10994cf0b 100644 --- a/starfish/test/test_overlap_utils.py +++ b/starfish/test/test_overlap_utils.py @@ -149,5 +149,5 @@ def test_take_max(): # The overlap section hits half of the spots from each intensity table, 5 from it1 # and 10 from i21. It2 wins and the resulting concatenated table should have all the - # spots from it2 (20) and 6 (one on the boarder) from it1 (6) for a total of 26 spots + # spots from it2 (20) and 6 (one on the border) from it1 (6) for a total of 26 spots assert concatenated.sizes[Features.AXIS] == 26 diff --git a/starfish/types/__init__.py b/starfish/types/__init__.py index ac7e55135..7ffb288f8 100644 --- a/starfish/types/__init__.py +++ b/starfish/types/__init__.py @@ -7,6 +7,7 @@ CORE_DEPENDENCIES, Features, LOG, + OverlapStrategy, PHYSICAL_COORDINATE_DIMENSION, PhysicalCoordinateTypes, STARFISH_EXTRAS_KEY diff --git a/starfish/util/overlap_utils.py b/starfish/util/overlap_utils.py index f4e423ff9..d6cb72e49 100644 --- a/starfish/util/overlap_utils.py +++ b/starfish/util/overlap_utils.py @@ -1,10 +1,10 @@ import itertools -from typing import List, Optional, Tuple +from typing import List, Optional, Sequence, Tuple import xarray as xr from starfish.types import Coordinates, Features, Number -from starfish.types._constants import OverlapStrategy +from starfish.types import OverlapStrategy class Area: @@ -18,7 +18,7 @@ def __init__(self, min_x: Number, max_x: Number, min_y: Number, max_y: Number): self.min_y = min_y self.max_y = max_y - def __eq__(self, other): + def __eq__(self, other) -> bool: return (self.min_x == other.min_x and self.min_y == other.min_y and self.max_x == other.max_x @@ -47,7 +47,7 @@ def find_intersection(area1: "Area", area2: "Area") -> Optional["Area"]: return None -def find_overlaps_of_xarrays(xarrays: List[xr.DataArray]) -> List[Tuple[int, int]]: +def find_overlaps_of_xarrays(xarrays: Sequence[xr.DataArray]) -> Sequence[Tuple[int, int]]: """ Find all the overlap areas within a list of xarrays. @@ -58,7 +58,7 @@ def find_overlaps_of_xarrays(xarrays: List[xr.DataArray]) -> List[Tuple[int, int Returns ------- - list[Tuple[int, int]] : + List[Tuple[int, int]] : A list of tuples containing the indices of two overlapping IntensityTables.