[RFC] Add linear unmixing

[kevinyamauchi]

Objective

This PR introduces a component to perform linear unmixing across channels to compensate for crosstalk between filters.

Overview

Per our discussion in #945, I implemented the following:

2. Channel cross-talk: takes a square channel x channel coefficient matrix. This uses the Filter PipelineComponent.

Usage

import numpy as np
from starfish.imagestack.imagestack import ImageStack

import starfish

# Create the ImageStack
im = np.ones((1, 3, 5, 2, 2))
stack = ImageStack.from_numpy_array(im)

# Create the coefficient matrix
coeff_mat = np.array([[1, 0, 0], [-0.25, 1, -0.25], [0, 0, 1]])

# Unmix
filter_unmix = starfish.image.Filter.LinearUnmixing(coeff_mat=coeff_mat)
stack2 = filter_unmix.run(stack, in_place=False, verbose=True)

Questions

• I think it would be nice to expose a couple of rescaling options. Perhaps normalizing to the max and/or mean?
• This assumes the unmixing will be performed across channels. Is there value in making it general? At the moment, I can't think of any use cases, so I'm leaning towards keeping it specific to channels.

[codecov-io]

Codecov Report

Merging #987 into master will increase coverage by 0.51%.
The diff coverage is 97.91%.

@@            Coverage Diff             @@
##           master     #987      +/-   ##
==========================================
+ Coverage   88.59%   89.11%   +0.51%     
==========================================
  Files         162      168       +6     
  Lines        5824     6514     +690     
==========================================
+ Hits         5160     5805     +645     
- Misses        664      709      +45

Impacted Files	Coverage Δ
starfish/test/image/filter/test_linear_unmixing.py	100% <100%> (ø)
starfish/image/_filter/linear_unmixing.py	96.96% <96.96%> (ø)
starfish/util/exec.py	85.36% <0%> (-5.26%)	⬇️
starfish/display/stack.py	23.61% <0%> (-3.81%)	⬇️
starfish/spots/_detector/_base.py	76.92% <0%> (-3.08%)	⬇️
starfish/spots/_detector/__init__.py	97.61% <0%> (-2.39%)	⬇️
...h/spots/_detector/trackpy_local_max_peak_finder.py	89.33% <0%> (-1.74%)	⬇️
...tarfish/test/spots/detector/test_spot_detection.py	100% <0%> (ø)	⬆️
starfish/starfish.py	100% <0%> (ø)	⬆️
...cent_features/test_intensities_to_decoded_image.py	100% <0%> (ø)	⬆️
... and 21 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 37c53d2...a07eca5. Read the comment docs.

[kevinyamauchi]

In terms of checking the functionality, I think the following test could be enough to get us started. I haven't used pytest. Is using an assert statements sufficient or are there some pytest-specific functions we should use? Other thoughts?

Also, would this go in starfish/test/pipeline?

import numpy as np
from starfish.imagestack.imagestack import ImageStack

import starfish

# Create image
im = np.ones((2, 3, 5, 2, 2))
stack = ImageStack.from_numpy_array(im)

# Create reference result
ref_result = np.ones((2, 3, 5, 2, 2))
ref_result[:, 1, ...] = 0.5 * np.ones((5, 2, 2))

# Unmix
coeff_mat = np.array([[1, 0, 0], [-0.25, 1, -0.25], [0, 0, 1]])
filter_unmix = starfish.image.Filter.LinearUnmixing(coeff_mat=coeff_mat)
stack2 = filter_unmix.run(stack, in_place=False, verbose=False)

# Check result
assert np.all(ref_result == stack2.xarray.values)

[ambrosejcarr]

This is great. I've left some comments, and responses to your questions:

• I think it would be nice to expose a couple of rescaling options. Perhaps normalizing to the max and/or mean?

We have some separate normalization components (scale by percentile, I'm going to add histogram matching soon). Mean could work. Max seems like it would fit to outliers? I think I'd keep these separate because there are good reasons for these to be more general than per-channel (see below).

• This assumes the unmixing will be performed across channels. Is there value in making it general? At the moment, I can't think of any use cases, so I'm leaning towards keeping it specific to channels.

I couldn't think of anything either. Keeping it channel specific makes sense to me.

Also, would this go in starfish/test/pipeline?

I'd put it in starfish/test/image/filter

[ambrosejcarr]

I think this would be great if it took an xarray with labeled axes. Then later you could make this change:

-        n_channels = image.shape[0]
-        for c in range(n_channels):
+        for c in image.coords[Axes.CH]:

[kevinyamauchi]

I agree that it would be nice to use the axis labels for safety, but the apply() passes a numpy array. I propose that we leave it as is for now and make and issue proposing that we change the apply() method to pass an xarray. Thoughts?

[ambrosejcarr]

Fixed in #1035

[ambrosejcarr]

The framing of B = AX makes me think there's a one-liner that captures the logic in this for loop. Does that seem crazy? I'd be happy to talk through this on a call. 👍

[kevinyamauchi]

My intuition was the same. However, I couldn't figure out how to do it without the for loop. Within the for loop (i.e., for a given channel), we could do something like:

np.dot(np.transpose(image[c, :, :]),  coeff_mat[c])

but originally, I thought my previous implementation was more clear. I'm definitely not a broadcasting wizard, so if you can figure out how to reshape the matrices to do it all in one line, I'd love to learn.

[ambrosejcarr]

In terms of checking the functionality, I think the following test could be enough to get us started. I haven't used pytest. Is using an assert statements sufficient or are there some pytest-specific functions we should use? Other thoughts?

pytest just uses assert. You're doing the right thing 👍

[kevinyamauchi]

@ambrosejcarr I've added a test and addressed your comments. Would you mind taking a look?

[ambrosejcarr]

Closed by #1045