Merge branch 'master' of https://github.com/NEUBIAS/training-resources

_includes/filter_neighbourhood/filter_variance_skimage_napari.py

@@ -1,46 +1,78 @@
-# %% [markdown]
-# ## Applying variance filters to an image
+# %% 
+# Apply variance filters to an EM image to aid foreground background segmentation
 
 # %%
-# Instantiate the napari viewer
+# Import libraries and
+# instantiate the napari viewer
+import numpy as np
 import napari
+import matplotlib.pyplot as plt
 from OpenIJTIFF import open_ij_tiff
 viewer = napari.Viewer()
 
 # %%
-# Read the intensity image and add it to napari
-image, axes, scales, units = open_ij_tiff('https://github.com/NEUBIAS/training-resources/raw/master/image_data/xy_8bit__em_mitochondria_er.tif')
+# Read the image and add it to napari
+# - Try to find an intensity threshold that separates the intracellular region from the bright background
+image, *_ = open_ij_tiff('https://github.com/NEUBIAS/training-resources/raw/master/image_data/xy_8bit__em_mitochondria_er.tif')
 viewer.add_image(image)
 
+# %% 
+# Plot an intensity histogram as another way to find a threshold
+plt.hist(image.flatten(), bins=np.arange(image.min(), image.max() + 1))
+
 # %%
-# Appreciate that one cannot readily segment the sample by a simple intensity threshold
+# Appreciate that a simple threshold does not suffice to distinguish the cytoplasm from the bright background outside the cell
 binary_image = image < 230
 viewer.add_image(binary_image)
 
 # %% 
-# Apply a variance filter 
+# Apply a variance filter to separate the two regions
+# - The idea is to make use of the fact that the background had less local intensity variations than the cytoplasm
+# - Note that a mean filter might also work here (please try), but we want to learn something new
 
+# %%
 # Convert to uint16 to be able to accommodate 
-# the square of the unit8 image
-import numpy as np
+# the square of the unit8 image, which is needed to represent the variance
+# - Appreciate that such datatype issues are common applying any kind of "image math"
 print(image.dtype)
 uint16_image = image.astype(np.uint16)
+print(uint16_image.dtype)
+
+# %% 
+# View the uint16 image in napari and check whether the datatype conversion changed the pixel values
+# - Going from a lower bit-depth to a higher bit-depth should normally not change the values
+viewer.add_image(uint16_image)
 
-# Use ndimage instead of skimage,
-# because skimage.filters.rank.mean 
-# only works up to unit16, which here specifically
-# would have been sufficient. However if the input image 
-# would have been uint16 we would have needed unit32
-# for the variance.
-from scipy import ndimage
-sqr_mean_image = ndimage.uniform_filter(uint16_image, (11, 11))**2
-mean_sqr_mean = ndimage.uniform_filter(uint16_image**2, (11, 11))
-var_image = mean_sqr_mean - sqr_mean_image
+# %%
+# We could not find a function to directly compute the local variance of an image,
+# thus we introduce a, very useful, generic approach to compute any kind of local neighborhood
+# filter
+
+# import the function
+from scipy.ndimage import generic_filter
+
+# define a function that computes one value from an array of data
+def local_variance(data):
+  value = np.var(data)
+  return value
+
+# apply this function locally to image
+var_image = generic_filter(uint16_image, function=local_variance, size=11)
+
+# %%
+# View the variance filtered the image to napari
+# Look for a threshold that separates the cell from the background
 viewer.add_image(var_image)
 
 # %%
-# Segment the foreground by applying a threshold
-# (Note that in this specific case a simple mean filter might have been sufficient
-# to achieve the same result, because the sample is on average darker than the background)
+# Also inspect the histogram to find a threshold
+# - It is interesting to compare the histograms of the variance images computed
+#   with different filter widths, to see how well there are two populations of intensity values
+plt.hist(var_image.flatten(), bins=np.arange(var_image.min(), var_image.max() + 1))
+plt.yscale('log')
+
+# %%
+# Segment the cellular region by applying a threshold to the variance filtered image
+# - It can be interesting to also try some auto-threshold method here
 binary_var_image = var_image > 100
 viewer.add_image(binary_var_image)

_includes/median_filter/median_filter_skimage_napari.py

@@ -1,52 +1,59 @@
-### Median Filtering
+# %%
+# Median filtering
 
+# %%
+# import modules
+import napari
 import numpy as np
+from skimage import filters
+from skimage.filters import rank
+from skimage.morphology import disk # Structuring element
+from OpenIJTIFF import open_ij_tiff
 
+# %%
 # Instantiate the napari viewer
-import napari
 viewer = napari.Viewer()
 
-# Read the intensity image
-#
-# (Replace the image path to explore the other example images)
-#
-from OpenIJTIFF import open_ij_tiff
-image1, axes1, scales1, units1 = open_ij_tiff('https://github.com/NEUBIAS/training-resources/raw/master/image_data/xy_8bit_binary__squares_different_size.tif')
-
-# Inspect image data type and values
-print('image type:', image1.dtype,'\n',
-      'image shape:', image1.shape,'\n',
-      'intensity min:',   np.min(image1),'\n',
-      'intensity max:',   np.max(image1),'\n'
-      )
-
-# View the intensity image
-viewer.add_image(image1, name='original image1')
-
-# First method using the median filter
-from skimage import filters
-# Second method using the mean filter
-from skimage.filters import rank
+# %%
+# Read and view the image
+image, *_ = open_ij_tiff('https://github.com/NEUBIAS/training-resources/raw/master/image_data/xy_8bit__two_noisy_squares_different_size.tif')
+viewer.add_image(image)
+
+# %% 
+# Compare small median and mean filter
+# - Appreciate that the median filter better preserves the edges
+median_1 = filters.median(image, disk(1))
+viewer.add_image(median_1)
+mean_1 = rank.mean(image, disk(1))
+viewer.add_image(mean_1)
+
+# %% 
+# Compare a larger median and mean filter
+# - Appreciate that a median filter eliminates the small square entirely
+median_3 = filters.median(image, disk(3))
+viewer.add_image(median_3)
+mean_3 = rank.mean(image, disk(3))
+viewer.add_image(mean_3)
+
+
+############################  New image  ################################
+
+# %%
+# Clear the napari viewer
+viewer.layers.clear()
+
+# %%
+# Load another image
+image, *_ = open_ij_tiff('https://github.com/NEUBIAS/training-resources/raw/master/image_data/xy_8bit__PCNA.tif')
+
+# %%
+# View the image and find the radius of the intra-nuclear speckles
+viewer.add_image(image)
+
+# %%
+# Remove small intra-nuclear structures 
+# - Appreciate that a median filter removes the intra-nuclear speckles while keeping the nuclear shape intact
+# - This can be helpful in many ways, e.g. some deep learning nuclear segmentation tools get confused by intra-nuclear speckles
+image_without_speckles = filters.median(image, disk(radius=15))
+viewer.add_image(image_without_speckles)
 
-from skimage.morphology import disk
-
-# Local median filtering with radius 1
-median_1 = filters.median(image1, disk(1))
-viewer.add_image(median_1, name='Median1')
-# Local mean filtering with radius 1
-mean_1 = rank.mean(image1, disk(1))
-viewer.add_image(mean_1, name='Mean1')
-
-# Local median filtering with radius 2
-Median_2 = filters.median(image1, disk(2))
-viewer.add_image(Median_2, name='Median2')
-# Local mean filtering with radius 2
-mean_2 = rank.mean(image1, disk(2))
-viewer.add_image(mean_2, name='Mean2')
-
-# Local median filtering with radius 5
-Median_3 = filters.median(image1, disk(5))
-viewer.add_image(Median_3, name='Median3')
-# Local mean filtering with radius 5
-mean_3 = rank.mean(image1, disk(5))
-viewer.add_image(mean_3, name='Mean3')