This page was generated from: notebooks/nanopores/how_to_segment_image_data_using_local_thresholds.ipynb
[1]:
%load_ext autoreload
%autoreload 2
%config InlineBackend.rc = {'figure.figsize': (10,6)}
%matplotlib inline
Segment image data using local thresholdsο
This notebook meshes an image with multiple domains using Nanomesh. The image contains a gradient from left to right. If the image background is relatively uniform, global thresholds can be used to separate the different domains. In this case, the image has a contrast gradient from left to right, so this example shows three different methods to deal with this.
The image is then meshed using Nanomesh to obtain a mesh with multiple domains. The mesh contains triangles labeled as βporeβ or βbulkβ material.
Loading and pre-processing the dataο
This example uses nanopore sample data from nanomesh.data
. Notice how the image has a bit of a gradient going from top to bottom. We apply a gaussian filter to reduce image noise.
If you want to use your own data, any numpy array can be passed to into a `Image
<https://nanomesh.readthedocs.io/en/latest/nanomesh.volume.html#nanomesh.volume.Volume>`__ object. Data stored as .npy
can be loaded using Image.load()
.
[2]:
from nanomesh import Image
from nanomesh.data import nanopores_gradient
data = nanopores_gradient()
plane = Image(data).gaussian(sigma=5)
plane.show()
[2]:
<AxesSubplot:xlabel='x', ylabel='y'>
Using .try_all_thresholds
is usually a good way to find a useful value to segment the data. In this case, the gradient prevents us from getting a useful segmentation.
[3]:
plane.try_all_threshold(figsize=(4, 10))
Local thresholdingο
This section explores three different methods for local thresholding.
User-defined local thresholdο
In a local threshold filter, the local neighbourhoud is used to define a threshold map for the image. The block_size
defines the size of the neighbourhood, a 101 by 101 pixel window in this case. The offset is used to tune the map. The advantage of this method is that it is fairly simple. The downside is that it requires manually tuning two parameters (offset
, blocksize
) to get a useful result.
[4]:
import matplotlib.pyplot as plt
offset = 150
block_size = 101
local_thresh = plane.threshold('local', block_size=block_size, offset=offset)
seg_local = plane > local_thresh
## Plot
import matplotlib.pyplot as plt
fig, axes = plt.subplots(ncols=2, figsize=(6, 6))
ax = axes.ravel()
ax[0].imshow(plane.image)
ax[0].set_title('Original')
ax[1].imshow(seg_local.image)
ax[1].set_title('Segmentation')
for a in ax:
a.axis('off')
plt.show()
Otsu local thresholdο
Applies an otsu rank filter to determine the local otsu threshold to segment on. Similar to the global otsu method, but constrained to a small area around each pixel. The advantage of this method that it only requires tuning a single parameter, namely the radius
of the neighbourhood. The downside is that only the otsu method is available as a rank filter.
[5]:
from skimage.morphology import disk
from skimage.filters import rank
from skimage.util import img_as_ubyte
plane_norm = plane.normalize_values().apply(img_as_ubyte)
radius = 41
selem = disk(radius)
local_otsu = plane_norm.apply(rank.otsu, selem=selem)
seg_local_otsu = plane_norm >= local_otsu
## Plot
import matplotlib.pyplot as plt
fig, axes = plt.subplots(ncols=3, figsize=(9, 6))
axes = axes.ravel()
axes[0].imshow(plane.image)
axes[0].set_title('Original')
axes[1].imshow(local_otsu.image)
axes[1].set_title(f'Local Otsu filter\n(radius={radius})')
axes[2].imshow(seg_local_otsu.image)
axes[2].set_title('Segmentation')
for ax in axes:
ax.axis('off')
plt.show()
Adaptive histogram equalizationο
This method tries to correct the image by removing the local gradient first, using adaptive histogram equalization. The advantage is that all global threshold finders are available. Another advantage is that there are no parameters to tune. At the same time, this is also a disadvantage if the result is not good π
[6]:
from skimage import exposure
plane_eq = plane.normalize_values()
plane_eq = plane_eq.apply(exposure.equalize_adapthist)
plane_eq.try_all_threshold(figsize=(6, 15))
Compare resultsο
The next cell creates a plot that compares the result of all three methods. For the histogram equalization, the Otsu and Li filters are shown.
[7]:
import matplotlib.pyplot as plt
fig, axes = plt.subplots(ncols=5, figsize=(9, 5))
axes = axes.ravel()
axes[0].imshow(plane.image)
axes[0].set_title('Original')
axes[1].imshow(seg_local.image)
axes[1].set_title('Local\nthresholding')
axes[2].imshow(seg_local_otsu.image)
axes[2].set_title('Local Otsu')
seg_clahe_otsu = plane_eq.binary_digitize(threshold='otsu')
axes[3].imshow(seg_clahe_otsu.image)
axes[3].set_title('Hist. equal.\n + global Otsu')
seg_clahe_li = plane_eq.binary_digitize(threshold='li')
axes[4].imshow(seg_clahe_li.image)
axes[4].set_title('Hist. equal.\n + global Li')
for ax in axes:
ax.axis('off')
plt.show()
Meshing the image dataο
Finally, itβs time to mesh the image.
[8]:
# Triangulation has some issues with small elements near the border,
# so these will be cleared first.
# https://github.com/hpgem/nanomesh/issues/86
seg_local = seg_local.clear_border(object_label=0, fill_val=1)
seg_local.show()
[8]:
<AxesSubplot:xlabel='x', ylabel='y'>
[9]:
%%time
from nanomesh import Mesher2D
mesher = Mesher2D(seg_local)
mesher.generate_contour(max_edge_dist=4)
mesh = mesher.triangulate(opts='q30a100')
Wall time: 164 ms
View the result using matplotlib
:
[10]:
mesh.plot(color_map={0: 'lightgray'}, lw=1)
[10]:
(<AxesSubplot:title={'center':'line mesh'}>,
<AxesSubplot:title={'center':'triangle mesh'}>)
Or, view the result using pyvista:
(Note: Use mesh.plot_itk()
for an interactive view).
[11]:
mesh.plot_pyvista(jupyter_backend='static', show_edges=False)
Save the data:
[12]:
mesh.write("mesh_gradient.msh", file_format='gmsh22', binary=False)
Warning: Appending zeros to replace the missing geometrical tag data.