Skip to content

Commit

Permalink
Update spot_detection
Browse files Browse the repository at this point in the history
  • Loading branch information
@psobolewskiPhD
psobolewskiPhD committed Apr 22, 2024
1 parent dc7ef57 commit c083608
Showing 1 changed file with 85 additions and 136 deletions.
221 changes: 85 additions & 136 deletions napari-workshops/notebooks/spot_detection.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ kernelspec:
name: python3
---

# Interactive analysis with Jupyter notebook, napari, scikit-image, and scipy
# Exploratory analysis: spot detection

## Overview
In this activity, we will perform spot detection on some in situ sequencing data
Expand All @@ -28,16 +28,23 @@ The data were downloaded from the
## Next steps

Following this activity, we will use the workflow generated in this activity to
create a napari spot detection plugin.
create a napari spot detection widget or even plugin.

## screenshots
For the solution notebook, we are including screenshots via the `nbscreenshot`
utility. These are not required for your notebook.
## `binder` setup
```{code-cell} ipython3
:tags: [remove-output]
# this cell is required to run these notebooks on Binder. Make sure that you also have a desktop tab open.
import os
if 'BINDER_SERVICE_HOST' in os.environ:
os.environ['DISPLAY'] = ':1.0'
```

An example usage: `nbscreenshot(viewer)`
## Screenshots
As previously, we will use the `nbscreenshot` to document our work in the notebook as we go along.
As a reminder, the usage is:

```{code-cell} python
from napari.utils import nbscreenshot
```Python
nbscreenshot(viewer)
```

## Load the data
Expand All @@ -60,42 +67,40 @@ spots = io.imread(spots_url)

## View the data

We will use napari to view our data. To do so, we first must create the viewer.
Once the Viewer is created, we can add images to the viewer via the Viewer's
`add_image()` method.

```{code-cell} ipython3
:tags: [remove-output]
# this cell is required to run these notebooks on Binder. Make sure that you also have a desktop tab open.
import os
if 'BINDER_SERVICE_HOST' in os.environ:
os.environ['DISPLAY'] = ':1.0'
```
To view our data, we will create a napari viewer and then we will add the images to the viewer
via the viewer's `add_image()` method. We will also import the `nbscreenshot` utility.

```{code-cell} python
import napari
from napari.utils import nbscreenshot
# create the napari viewer
viewer = napari.Viewer();
viewer = napari.Viewer()
# add the nuclei image to the viewer
viewer.add_image(nuclei);
viewer.add_image(nuclei, colormap='green')
# add the spots image to the viewer
viewer.add_image(spots, colormap='magenta', blending='additive')
```

In the cell below, add the spots image to the viewer as was done above for the
nuclei image. After loading the data, inspect it in the viewer and adjust the
layer settings to your liking (e.g., contrast limits, colormap). You can
pan/zoom around the image by click/dragging to pan and scrolling with your
After loading the data, inspect it in the viewer and adjust the
layer settings to your liking (e.g., contrast limits, colormap).
You can pan/zoom around the image by click/dragging to pan and scrolling with your
mousewheel or trackpad to zoom.

```{tip}
You can adjust a layer's opacity to see the change how much you see of the
layers that are "under" it.
```


Once you are satisfied, you can print the output of any manual changes and then take a screenshot of the viewer.

```{code-cell} python
# add the spots image to the viewer
viewer.add_image(spots)
print('Colormap: ', viewer.layers['nuclei'].colormap)
print('Contrast limits: ', viewer.layers['nuclei'].contrast_limits)
print('Opacity: ', viewer.layers['nuclei'].opacity)
```

```{code-cell} python
Expand All @@ -104,141 +109,84 @@ nbscreenshot(viewer)

## Create an image filter

You may have noticed the the spots image contains background and
autofluorescence from the cells. To improve spot detection, we will apply a high
pass filter to improve the contrast of the spots.
If you look carefull at the `spots` layer, you will notice that it contains background and
autofluorescence from the cells. It may help to just look at the single channel and for better
contrast you could choose an inverted colormap (and minimum blending if you want to view
multiple layers).

```{code-cell} python
viewer.layers['nuclei'].visible = False
viewer.layers['spots'].colormap = "I Orange"
```

To improve spot detection, we will apply a high pass filter to improve the contrast of the spots.

A simple way to achieve this is to:
1. blur the image with a Gaussian, which removes high-frequency information—this is why we use it for denoising.
2. subtract the blurred image from the original.

Lets try this using the `gaussian_filter` from `scipy` with a sigma of 2 px and lets clip any negative values:

```{code-cell} python
import numpy as np
from scipy import ndimage as ndi
def gaussian_high_pass(image: np.ndarray, sigma: float = 2):
"""Apply a gaussian high pass filter to an image.
Parameters
----------
image : np.ndarray
The image to be filtered.
sigma : float
The sigma (width) of the gaussian filter to be applied.
The default value is 2.
Returns
-------
high_passed_im : np.ndarray
The image with the high pass filter applied
"""
low_pass = ndi.gaussian_filter(image, sigma)
high_passed_im = image - low_pass
return high_passed_im
low_pass = ndi.gaussian_filter(spots, 2)
high_passed_spots = (spots - low_pass).clip(0)
```

In the cell below, apply the gaussian high pass filter to the `spots` image and
add the image to the viewer.
Let's visualize both versions of the data.

```{code-cell} python
# Use the gaussian_high_pass function to filter the spots image
filtered_spots = gaussian_high_pass(spots, 2)
# add the filtered image to the viewer
# hint: set the opacity < 1 in order to see the layers underneath
viewer.add_image(filtered_spots, opacity=0.6, colormap='viridis')
viewer.add_image(low_pass, colormap="I Forest", blending="minimum")
viewer.add_image(high_passed_spots, colormap="I Blue", blending="minimum")
```

Toggle the visibility of the 3 spots layers to get a better sense of the effect of our filter.
Let's hide the `low_pass` and take a screenshot for documentation.


```{code-cell} python
viewer.layers['low_pass'].visibility = False
nbscreenshot(viewer)
```

## Detect spots

Next, we will create a function to detect the spots in the spot image. This
function should take the raw image, apply the gaussian high pass filter from
above and then use one of the blob detection algorithms from sci-kit image to
perform the blob detection. The `detect_spots()` function should return a numpy
array containing the coordinates of each spot and a numpy array containing the
diameter of each spot.
Next, we will use one of the blob detection algorithms from sci-kit image to
perform the blob detection.

Some hints:
- See the [blob detection tutorial from scikit-image](https://scikit-image.org/docs/dev/auto_examples/features_detection/plot_blob.html). We recommend the [blob_log detector](https://scikit-image.org/docs/dev/api/skimage.feature.html#skimage.feature.blob_log), but feel free to experiment!
```{tip}
- See the [blob detection tutorial from scikit-image](https://scikit-image.org/docs/stable/auto_examples/features_detection/plot_blob.html). We recommend the [blob_log detector](https://scikit-image.org/docs/stable/api/skimage.feature.html#skimage.feature.blob_log), but feel free to experiment!
- See the "Note" from the blob_log docs: "The radius of each blob is approximately $\sqrt{2}\sigma$ for a 2-D image"
```

```{code-cell} python
import numpy as np
from skimage.feature import blob_log
def detect_spots(
image: np.ndarray,
high_pass_sigma: float = 2,
spot_threshold: float = 0.01,
blob_sigma: float = 2
):
"""Apply a gaussian high pass filter to an image.
Parameters
----------
image : np.ndarray
The image in which to detect the spots.
high_pass_sigma : float
The sigma (width) of the gaussian filter to be applied.
The default value is 2.
spot_threshold : float
The threshold to be passed to the blob detector.
The default value is 0.01.
blob_sigma: float
The expected sigma (width) of the spots. This parameter
is passed to the "max_sigma" parameter of the blob
detector.
Returns
-------
points_coords : np.ndarray
An NxD array with the coordinate for each detected spot.
N is the number of spots and D is the number of dimensions.
sizes : np.ndarray
An array of size N, where N is the number of detected spots
with the diameter of each spot.
"""
# filter the image
filtered_spots = gaussian_high_pass(image, high_pass_sigma)
# detect the spots on the filtered image
blobs_log = blob_log(
filtered_spots,
max_sigma=blob_sigma,
num_sigma=1,
threshold=spot_threshold
)
# detect the spots on the filtered image
blobs_log = blob_log(
high_passed_im,
max_sigma=3,
threshold=None, # use a relative threshold instead
threshold_rel=0.2
)
# convert the output of the blob detector to the
# desired points_coords and sizes arrays
# (see the docstring for details)
points_coords = blobs_log[:, 0:2]
sizes = 2 * np.sqrt(2) * blobs_log[:, 2]
return points_coords, sizes
# convert the output of the blob detector to the
# desired points_coords and sizes arrays
# (see the docstring for details)
spot_coords = blobs_log[:, 0:2]
spot_sizes = 2 * np.sqrt(2) * blobs_log[:, 2]
```

In the cell below, apply `detect_spots()` to our `spots` image. To visualize the
results, add the spots to the viewer as a
[Points layer](https://napari.org/tutorials/fundamentals/points.html). If you
To visualize the results, add the spots to the viewer as a
[Points layer](https://napari.org/stable/tutorials/fundamentals/points.html). If you
would like to see an example of using a points layer, see
[this example](https://napari.org/gallery/add_points.html). To test out your
function, vary the detection parameters and see how they affect the results.
Note that each time you run the cell, the new results are added as an addition
Points layer, allowing you to compare results from different parameters.
[this example](https://napari.org/gallery/add_points.html).

```{code-cell} python
# detect the spots
spot_coords, spot_sizes = detect_spots(
spots,
high_pass_sigma=2,
spot_threshold=0.01,
blob_sigma=2
)
# add the detected spots to the viewer as a Points layer
viewer.add_points(spot_coords, size=spot_sizes)
```
Expand All @@ -248,7 +196,8 @@ nbscreenshot(viewer)
```

## Conclusion
In this activity, we have interactively prototyped a spot detection function
using a combination of jupyter notebook, scipy, scikit-image, and napari. In the
next activity, we will take the spot detection function we created and turn it
into a napari plugin.
In this activity, we have developed done the exploratory analysis for spot detection
using a combination of Jupyter notebook, scipy, scikit-image, and napari. In the
next activity, we will convert this workflow into a spot detection function to make
it easier to explore the effect of parameters. Then, we will turn that into a napari widget
using `magicgui`.

0 comments on commit c083608

Please sign in to comment.