Mask Pixels when carrying out pixel detection and decoding

[ambrosejcarr]

The Expansion Microscopy Pipeline uses a filter to mask pixels that are ineligible for spot finding. This filter:

1. Matches the intensity distribution of each channel in each round
2. Averages the channels in each round, producing a single gray-scale volume per round
3. Averages the rounds, producing a single gray-scale volume per experiment
4. Thresholds this image to identify eligible pixels for decoding

We need to implement this filter to support expansion microscopy.

Note: because this reduces the overall pixels eligible for decoding, it may alleviate memory/cpu requirements in #1282

[kevinyamauchi]

@shanaxel42 patiently discussed this with me yesterday and based on that, I think the best approach here is to use the segmentation mask to select the pixels corresponding to spots from the IntensityTable. Rather than making a monolithic filter component, I think we can make/modify a few components, which will give users more flexibility in the future.

Proposed pipeline

Assumptions

1. We have an ImageStack, im_stack that is fully registered and filtered as desired. Each channel corresponds to a nucleotide.
2. We have a segmentation component that can operate on ImageStacks without a nuclei image. The resulting segmentation mask is an xarray with coordinates/dims that match im_stack

Pipeline

# Sum along the round and channel axes. (To do  2)
proj_stack = im_stack.project('sum', ('Axes.ROUND, Axes.CH'))

# Segment (To do 3)
seg = Segment.Watershed()
masks = seg.run(proj_stack)

# Make a single label image (xarray) from the SegmentationMaskCollection (to do 4)
single_mask = mask.single_image()

# Make the IntensityTable
it = IntensityTable.from_image_stack(im_stack)

# Use the segmentation mask to select the pixels corresponding to spots
it_masked = it.where(single_mask, drop=True)

# Decode it_masked with a pixel-wise decoder

Potential variants

• Masks can also be generated via spot finding. We could make a component that generates a mask from a set of spot coordinates and radii.
• We could potentially perform the decoding on a per spot basis (loop through the SegmentedMaskCollection

To do

1. Verify the projected mask (dims=(x, y) or (x, y, z)) broadcasts to (r, c, z, x, y) such it_masked = it.where(single_mask, drop=True) works.
2. Create a method to project along an axis with averaging or summing (#1302)
3. Either create a new segmentation component or modify Segment.Watershed such that it can accept other things for seeds than an image of nuclei.
4. Add method to output a single mask image from SegmentedMaskCollection

Thoughts, @ambrosejcarr, @dganguli ?

[ambrosejcarr]

I like this. One concern I have is how this will function for large fields of view that are dense in z, for example expansion microscopy, whose fields of view are (4, 4, 300, 2048, 2048).

For these images, it would be great to be able to mask individual volumes (per round) and combine the results, as this would substantially reduce memory usage and make processing more parallel.

I outline my best guess at this workflow in #876.

How would the solution proposed here interact with the goals I suggest in #876?

[kevinyamauchi]

I agree that that the process described above wouldn't work well with acquisitions densely sampled in Z. However, I think the general approach is totally compatible by just performing the operations on a per-round basis. In terms of the operations, it is the same as if it were done in a monolithic pipeline component. However, it would require the user to write a bit more code (e.g.,, some sort of loop or multiprocessing map).

For example, if you are okay with doing the spot segmentation on a per-round basis (I think this is the same as the approach you described in #876):

1. perform the projection and segmentation I proposed above round-by-round (i.e., select one round from the full field of view and proceed)
2. concatenate the masked intensity tables
3. decode.

Alternatively, if you'd like to perform the segmentation on all rounds:

1. perform the (round, channel) projection on the full FOV, which generates a much smaller image
2. generate the spot mask
3. mask on a per-round basis with the spot mask generated in (2)
4. concatenate the masked intensity tables
5. decode.

[ambrosejcarr]

I think your proposed implementation makes sense. This looks good to implement to me. Making sure I'm understanding, this will entail:

1. A segmentation component
2. A modification of the PixelSpotDecoder that allows it to take as input a mask from (1) above

For 1 & 3 to work, I think the base calling filter should optionally take a codebook and translate it to the new base-called feature space. I don't think we did that work yet, right?

[kevinyamauchi]

Hey @ambrosejcarr ! Sorry it's taken so long. I ended up needing to dig into the PixelSpotDetector more than I thought. In short, just masking the ImageStack or IntensityTable and feeding it into the PixelSpotDetector is not straightforward because CombineAdjacentFeatures reconstructs an image by reshaping the decoded IntensityTable. Thus, non-rectangular masks don't work. Also, the pixel coordinates for the spots are determined from the reshaped image (a numpy array) and not the IntensityTable, so the pixel coordinates of the masked features are not preserved.

I propose we do one of the following:

1. Replace masked out pixels with 0 and run the decoding on all pixels. This is a straightforward change, but I don't think really achieves the spirit of this issue (decoding is still caried out on all pixels)
2. Drop the masked out pixels from the IntensityTable before decoding, but place them back into the full FOV volume in the decoded_image. This would eliminate decoding of masked out pixels, but the connected components analysis would be carried out on the full image volume, so I'm not sure it would get the full speedup.
3. Convert the mask to a SegmentationMaskCollection and perform pixel wise decoding on each mask region. We would iterate through the segmentation mask collection and perform pixel-wise decoding on the selected pixels in each mask region, then concatenate the features into an IntensityTable. This would ensure compatibility with the starfish segementation components and allow users to generate their own masks.

It seems like (1) is only useful if you need something really fast just to get started. Options (2) and (3) seem pretty similar and I think option (3) probably depends on option (2). Therefore, I think we should do option (2) now and then add (3) in a later PR if needed. I am also open to other approaches. Thoughts?