Repository: spacetx/starfish

Issue:
`IntensityTable` is truncating x,y,z coordinates for spot detection methods that return subpixel results
We've been in the process of comparing the results of our current pipeline vs. a starfish version and it looks like we've stumbled on a really unfortunate 'feature' of the current `IntensityTable` data structure for spot based detection methods.

Below is a visual example (diamonds are our pipeline, circles are starfish), with the takeaway being that the starfish circles are not centered on detected spots:
![Screenshot from 2019-05-02 14-44-25](https://user-images.githubusercontent.com/43766899/57182459-2941d800-6e54-11e9-8439-17f0be3da602.png)


Both our spot finder and the [Starfish equivalent](http://soft-matter.github.io/trackpy/dev/tutorial/uncertainty.html) are based on the Crocker-Grier centroid finding algorithm and can determine the coordinates of a spot with sub-pixel resolution.

So I ended up chasing a [red herring](https://github.com/spacetx/starfish/pull/1300) before discovering that `IntensityTable` x, y, z properties are `int64`, which would `floor()` any result returned by the Crocker-Grier algorithm. In fact, this flooring appears to happen [immediately after results are produced](https://github.com/spacetx/starfish/blob/master/starfish/core/spots/_detect_spots/trackpy_local_max_peak_finder.py#L137).

It looks like (xc, yc, zc) are then real world coordinates that map to each set of pixel indices (x, y, z). I think that this design is probably (?) fine for pixel based decoding methods, but throwing out sub-pixel information for spot methods is a big deal breaker.

CC: @berl @ambrosejcarr 

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

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

```diff
@@            Coverage Diff             @@
##           master    #1330      +/-   ##
==========================================
+ Coverage   88.96%   88.98%   +0.01%     
==========================================
  Files         146      146              
  Lines        5258     5257       -1     
==========================================
  Hits         4678     4678              
+ Misses        580      579       -1
```


| [Impacted Files](https://codecov.io/gh/spacetx/starfish/pull/1330?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [...ots/\_detect\_spots/trackpy\_local\_max\_peak\_finder.py](https://codecov.io/gh/spacetx/starfish/pull/1330/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9zcG90cy9fZGV0ZWN0X3Nwb3RzL3RyYWNrcHlfbG9jYWxfbWF4X3BlYWtfZmluZGVyLnB5) | `91.22% <ø> (-0.44%)` | :arrow_down: |
| [starfish/core/spots/\_detect\_spots/blob.py](https://codecov.io/gh/spacetx/starfish/pull/1330/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9zcG90cy9fZGV0ZWN0X3Nwb3RzL2Jsb2IucHk=) | `96.15% <100%> (ø)` | :arrow_up: |
| [starfish/core/intensity\_table/intensity\_table.py](https://codecov.io/gh/spacetx/starfish/pull/1330/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9pbnRlbnNpdHlfdGFibGUvaW50ZW5zaXR5X3RhYmxlLnB5) | `96.42% <100%> (ø)` | :arrow_up: |
| [...ore/intensity\_table/intensity\_table\_coordinates.py](https://codecov.io/gh/spacetx/starfish/pull/1330/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS9pbnRlbnNpdHlfdGFibGUvaW50ZW5zaXR5X3RhYmxlX2Nvb3JkaW5hdGVzLnB5) | `100% <100%> (ø)` | :arrow_up: |
| [starfish/core/types/\_constants.py](https://codecov.io/gh/spacetx/starfish/pull/1330/diff?src=pr&el=tree#diff-c3RhcmZpc2gvY29yZS90eXBlcy9fY29uc3RhbnRzLnB5) | `100% <0%> (+1.88%)` | :arrow_up: |

------

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


Patch:
diff --git a/starfish/core/intensity_table/intensity_table.py b/starfish/core/intensity_table/intensity_table.py
index 30951497..8f6ab949 100644
--- a/starfish/core/intensity_table/intensity_table.py
+++ b/starfish/core/intensity_table/intensity_table.py
@@ -285,8 +285,14 @@ class IntensityTable(xr.DataArray):
 
     @classmethod
     def synthetic_intensities(
-            cls, codebook, num_z: int=12, height: int=50, width: int=40, n_spots=10,
-            mean_fluor_per_spot=200, mean_photons_per_fluor=50
+            cls,
+            codebook,
+            num_z: int = 12,
+            height: int = 50,
+            width: int = 40,
+            n_spots: int = 10,
+            mean_fluor_per_spot: int = 200,
+            mean_photons_per_fluor: int = 50,
     ) -> "IntensityTable":
         """
         Creates an IntensityTable with synthetic spots, that correspond to valid
@@ -317,8 +323,9 @@ class IntensityTable(xr.DataArray):
 
         # TODO nsofroniew: right now there is no jitter on x-y positions of the spots
         z = np.random.randint(0, num_z, size=n_spots)
-        y = np.random.randint(0, height, size=n_spots)
-        x = np.random.randint(0, width, size=n_spots)
+        y = np.random.uniform(0, height - 1, size=n_spots)
+        x = np.random.uniform(0, width - 1, size=n_spots)
+
         r = np.empty(n_spots)
         r.fill(np.nan)  # radius is a function of the point-spread gaussian size
         spot_attributes = SpotAttributes(
diff --git a/starfish/core/intensity_table/intensity_table_coordinates.py b/starfish/core/intensity_table/intensity_table_coordinates.py
index 87737ed3..2b27b1df 100644
--- a/starfish/core/intensity_table/intensity_table_coordinates.py
+++ b/starfish/core/intensity_table/intensity_table_coordinates.py
@@ -33,11 +33,15 @@ def transfer_physical_coords_from_imagestack_to_intensity_table(
         return intensity_table
 
     for axis, coord in pairs:
-
-        imagestack_pixels: np.ndarray = image_stack.xarray[axis].values
-        intensity_table_pixels: np.ndarray = intensity_table[axis].values.astype(int)
-        pixel_inds: np.ndarray = imagestack_pixels[intensity_table_pixels]
-        coordinates: xr.DataArray = image_stack.xarray[coord][pixel_inds]
+        imagestack_pixels: xr.DataArray = image_stack.xarray[axis]
+        intensity_table_pixel_offsets: np.ndarray = intensity_table[axis].values
+
+        # can't interpolate if the axis size == 1, so just select in that case.
+        if len(imagestack_pixels) == 1:
+            coordinate_fetcher = imagestack_pixels.sel
+        else:
+            coordinate_fetcher = imagestack_pixels.interp
+        coordinates = coordinate_fetcher({axis: intensity_table_pixel_offsets})[coord]
 
         intensity_table[coord] = xr.DataArray(coordinates.values, dims='features')
 
diff --git a/starfish/core/spots/_detect_spots/blob.py b/starfish/core/spots/_detect_spots/blob.py
index 4109fda8..14a771ed 100644
--- a/starfish/core/spots/_detect_spots/blob.py
+++ b/starfish/core/spots/_detect_spots/blob.py
@@ -115,7 +115,7 @@ class BlobDetector(DetectSpotsAlgorithmBase):
         fitted_blobs[Features.SPOT_RADIUS] = converted_radius
 
         # convert the array to int so it can be used to index
-        rounded_blobs = SpotAttributes(fitted_blobs.astype(int))
+        rounded_blobs = SpotAttributes(fitted_blobs)
 
         rounded_blobs.data['intensity'] = measure_spot_intensity(
             data_image, rounded_blobs, self.measurement_function)
diff --git a/starfish/core/spots/_detect_spots/trackpy_local_max_peak_finder.py b/starfish/core/spots/_detect_spots/trackpy_local_max_peak_finder.py
index be3a5354..6e614b36 100644
--- a/starfish/core/spots/_detect_spots/trackpy_local_max_peak_finder.py
+++ b/starfish/core/spots/_detect_spots/trackpy_local_max_peak_finder.py
@@ -133,10 +133,6 @@ class TrackpyLocalMaxPeakFinder(DetectSpotsAlgorithmBase):
             attributes.columns = new_colnames
 
         attributes['spot_id'] = np.arange(attributes.shape[0])
-        # convert these to int so it can be used to index
-        attributes.x = attributes.x.astype(int)
-        attributes.y = attributes.y.astype(int)
-        attributes.z = attributes.z.astype(int)
         return SpotAttributes(attributes)
 
     def run(


Test patch:
diff --git a/starfish/core/intensity_table/test/test_intensity_table_coords.py b/starfish/core/intensity_table/test/test_intensity_table_coords.py
index cf14eea1..be233735 100644
--- a/starfish/core/intensity_table/test/test_intensity_table_coords.py
+++ b/starfish/core/intensity_table/test/test_intensity_table_coords.py
@@ -1,3 +1,4 @@
+import math
 import random
 from collections import OrderedDict
 from typing import Union
@@ -57,13 +58,17 @@ def test_tranfering_physical_coords_to_intensity_table():
     # Assert that the physical coords align with their corresponding pixel coords
     for spot in xc.features:
         pixel_x = spot[Axes.X.value].data
-        physical_x = stack.xarray[Coordinates.X.value][pixel_x]
-        assert np.isclose(spot[Coordinates.X.value], physical_x)
+        pixel_x_floor, pixel_x_ceiling = math.floor(pixel_x), math.ceil(pixel_x)
+        physical_x_floor = stack.xarray[Coordinates.X.value][pixel_x_floor]
+        physical_x_ceiling = stack.xarray[Coordinates.X.value][pixel_x_ceiling]
+        assert physical_x_floor <= spot[Coordinates.X.value] <= physical_x_ceiling
 
     for spot in yc.features:
         pixel_y = spot[Axes.Y.value].data
-        physical_y = stack.xarray[Coordinates.Y.value][pixel_y]
-        assert np.isclose(spot[Coordinates.Y.value], physical_y)
+        pixel_y_floor, pixel_y_ceiling = math.floor(pixel_y), math.ceil(pixel_y)
+        physical_y_floor = stack.xarray[Coordinates.Y.value][pixel_y_floor]
+        physical_y_ceiling = stack.xarray[Coordinates.Y.value][pixel_y_ceiling]
+        assert physical_y_floor <= spot[Coordinates.Y.value] <= physical_y_ceiling
 
     # Assert that zc value is middle of z range
     for spot in zc.features:
diff --git a/starfish/core/intensity_table/test/test_synthetic_intensities.py b/starfish/core/intensity_table/test/test_synthetic_intensities.py
index 10c12160..6f35cb35 100644
--- a/starfish/core/intensity_table/test/test_synthetic_intensities.py
+++ b/starfish/core/intensity_table/test/test_synthetic_intensities.py
@@ -18,7 +18,7 @@ def test_synthetic_intensity_generation():
     a known random seed, that the output spots decode to match a target in the input Codebook
     """
     # set seed to check that codebook is matched. This seed generates 2 instances of GENE_B
-    np.random.seed(1)
+    np.random.seed(2)
     codebook = codebook_array_factory()
     num_z, height, width = 3, 4, 5
     intensities = IntensityTable.synthetic_intensities(
@@ -26,7 +26,7 @@ def test_synthetic_intensity_generation():
         num_z=num_z,
         height=height,
         width=width,
-        n_spots=2
+        n_spots=2,
     )
 
     # sizes should match codebook
diff --git a/starfish/core/test/factories.py b/starfish/core/test/factories.py
index b2d814f5..26f273f8 100644
--- a/starfish/core/test/factories.py
+++ b/starfish/core/test/factories.py
@@ -214,7 +214,7 @@ class SyntheticData:
                 warnings.simplefilter('ignore', FutureWarning)
                 values = spots.where(spots, drop=True)
 
-            image[round_, ch, values.z, values.y, values.x] = values
+            image[round_, ch, values.z, values.y.astype(int), values.x.astype(int)] = values
 
         intensities.values = img_as_float32(intensities)