Papers / Sites¶
Voronoi Tesselations
Ghosh, S. (1997). Tessellation-based computational methods for the characterization and analysis of heterogeneous microstructures. Composites Science and Technology, 57(9-10), 1187–1210
Self-Avoiding / Nearest Neighbor
Schwarz, H., & Exner, H. E. (1983). The characterization of the arrangement of feature centroids in planes and volumes. Journal of Microscopy, 129(2), 155–169.
Kubitscheck, U. et al. (1996). Single nuclear pores visualized by confocal microscopy and image processing. Biophysical Journal, 70(5), 2067–77.
Alignment / Distribution Tensor
Mader, K. et al (2013). A quantitative framework for the 3D characterization of the osteocyte lacunar system. Bone, 57(1), 142–154
Aubouy, M., et al. (2003). A texture tensor to quantify deformations. Granular Matter, 5, 67–70. Retrieved from http://arxiv.org/abs/cond-mat/0301018
Two point correlation
Dinis, L., et. al. (2007). Analysis of 3D solids using the natural neighbour radial point interpolation method. Computer Methods in Applied Mechanics and Engineering, 196(13-16)
In [1]:
import seaborn as sns
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (8, 8)
plt.rcParams["figure.dpi"] = 150
plt.rcParams["font.size"] = 14
plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = ['DejaVu Sans']
plt.style.use('ggplot')
sns.set_style("whitegrid", {'axes.grid': False})
Previously on QBI ...¶
- Image Enhancment
- Highlighting the contrast of interest in images
- Minimizing Noise
- Understanding image histograms
- Automatic Methods
- Component Labeling
- Single Shape Analysis
- Complicated Shapes (Thickness Maps)
Outline¶
- Motivation (Why and How?)
- Skeletons
- Tortuosity
- Watershed Segmentation
- Connected Objects
Metrics¶
We examine a number of different metrics in this lecture and additionally to classifying them as Local and Global we can define them as point and voxel-based operations.
Point Operations¶
- Nearest Neighbor
- Point (Center of Volume)-based Voronoi Tesselation
- Alignment
Voxel Operation¶
- Voronoi Tesselation
- Neighbor Counting
- 2-point (N-point) correlation function
In [2]:
from skimage.morphology import binary_opening, binary_closing, disk
from scipy.ndimage import distance_transform_edt
import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread
%matplotlib inline
bw_img = imread("ext-figures/bonegfiltslice.png")[::2, ::2]
thresh_img = binary_closing(binary_opening(bw_img < 90, disk(1)), disk(2))
fg_dmap = distance_transform_edt(thresh_img)
bg_dmap = distance_transform_edt(1-thresh_img)
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(20, 6), dpi=100)
ax1.imshow(bw_img, cmap='bone')
ax2.imshow(thresh_img, cmap='bone')
ax3.set_title('Segmentation')
ax3.imshow(fg_dmap, cmap='nipy_spectral')
ax3.set_title('Distance Map\nForeground')
ax4.imshow(bg_dmap, cmap='nipy_spectral')
ax4.set_title('Distance Map\nBackground')
Out[2]:
Text(0.5, 1.0, 'Distance Map\nBackground')
Distribution Objectives - finding the right questions¶
We want to know how many cells are alive¶
- Maybe small cells are dead and larger cells are alive $\rightarrow$ examine the volume distribution
- Maybe living cells are round and dead cells are really spiky and pointy $\rightarrow$ examine anisotropy
We want to know where the cells are alive or most densely packed¶
- We can visually inspect the sample (maybe even color by volume)
- We can examine the raw positions (x,y,z) but what does that really tell us?
- We can make boxes and count the cells inside each one
- How do we compare two regions in the same sample or even two samples?
So what do we still need¶
- A way for counting cells in a region and estimating density without creating arbitrary boxes
- A way for finding out how many cells are near a given cell, it's nearest neighbors
- A way for quantifying how far apart cells are and then comparing different regions within a sample
- A way for quantifying and comparing orientations
What would be really great?¶
A tool which could be adapted to answering a large variety of problems
- multiple types of structures
- multiple phases
Destructive Measurements¶
With most imaging techniques and sample types, the task of measurement itself impacts the sample.
- Even techniques like X-ray tomography which claim to be non-destructive still impart significant to lethal doses of X-ray radition for high resolution imaging
- Electron microscopy, auto-tome-based methods, histology are all markedly more destructive and make longitudinal studies impossible
- Even when such measurements are possible
- Registration can be a difficult task and introduce artifacts
Why is this important?¶
- techniques which allow us to compare different samples of the same type.
- are sensitive to common transformations
- Sample B after the treatment looks like Sample A stretched to be 2x larger
- The volume fraction at the center is higher than the edges but organization remains the same
Structure analysis¶
What we want to know about networks¶
- which structures are connected
- how they are connected
- Are there loops
- express the network in a simple manner
- quantify tortuosity
- branching
We start with a simpler example from the EPFL Dataset: EPFL CVLab's Library of Tree-Reconstruction Examples (http://cvlab.epfl.ch/data/delin)
An example - a street¶
In [60]:
import matplotlib.pyplot as plt # for showing plots
from skimage.io import imread # for reading images
import numpy as np # for matrix operations and array support
import pandas as pd # for reading the swc files (tables of somesort)
def read_swc(in_path):
swc_df = pd.read_csv(in_path, sep=' ', comment='#',
header=None)
# a pure guess here
swc_df.columns = ['id', 'junk1', 'x', 'y', 'junk2', 'width', 'next_idx']
return swc_df[['x', 'y', 'width']]
im_data = imread('ext-figures/ny_7.tif')
mk_data = read_swc('ext-figures/ny_7.swc')
fig, (ax1, ax3) = plt.subplots(1, 2, figsize=(15, 6))
ax1.imshow(im_data);ax1.set_title('Aerial Image')
ax3.imshow(im_data, cmap='bone') ;ax3.scatter(mk_data['x'], mk_data['y'], s=mk_data['width'], alpha=0.5,color='red',label="Roads") ;ax3.set_title('Annotated image'), ax3.legend();
A close-up of the street¶
In [61]:
im_crop = im_data[250:420:1, 170:280:1]
mk_crop = mk_data.query('y>250').query(
'y<420').query('x>170').query('x<280').copy()
mk_crop.x = (mk_crop.x-170)/1
mk_crop.y = (mk_crop.y-250)/1
fig, (ax1, ax3) = plt.subplots(1, 2, figsize=(15, 6))
ax1.imshow(im_crop)
ax1.set_title('Aerial Image')
ax3.imshow(im_crop, cmap='bone')
ax3.scatter(mk_crop['x'], mk_crop['y'], s=mk_crop['width'],color='red', alpha=0.25)
ax3.set_title('Roads')
Out[61]:
Text(0.5, 1.0, 'Roads')
In [62]:
import seaborn as sns
from skimage.morphology import opening, closing, disk # for removing small objects
from skimage.color import rgb2hsv
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16, 6))
def thresh_image(in_img):
v_img = rgb2hsv(in_img)[:, :, 2]
th_img = v_img > 0.4
op_img = opening(th_img, disk(1))
return op_img
ax1.imshow(im_crop)
ax1.set_title('Aerial Image')
ax2.imshow(rgb2hsv(im_crop)[:, :, 2],
cmap='bone')
ax2.set_title('HSV Representation')
seg_img = thresh_image(im_crop)
ax3.imshow(seg_img, cmap='bone')
ax3.set_title('Segmentation');
Other ways to do the segmentation...¶
Here we used a color space transformation
We could also use:
- Unsupervised segmentation like k-means
- Supervised segmentation like k-Nearest Neighbors (requires training)
Identify structures¶
Using connected component labeling
In [63]:
import seaborn as sns
from skimage.morphology import label
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 7))
ax1.imshow(im_crop)
ax1.set_title('Aerial Image')
lab_img = label(seg_img)
ax2.imshow(lab_img, cmap='nipy_spectral')
ax2.set_title('Labeling')
sns.heatmap(lab_img[::6, ::6], # Show every 6th pixel
annot=True,
fmt="d",
cmap='flag',
ax=ax3,
cbar=False,
vmin=0,
vmax=lab_img.max(),
annot_kws={"size": 8})
ax3.set_title('Labels');
Skeletonization - take one¶
The first step is to take the distance transform the structure
$$I_d(x,y) = \textrm{dist}(I(x,y))$$We can see in this image there are already local maxima that form a sort of backbone which closely maps to what we are interested in.
In [64]:
from scipy import ndimage
keep_lab_img = lab_img == 1
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 7))
ax1.imshow(keep_lab_img)
ax1.set_title('Road Component\n(Largest)')
dist_map = ndimage.distance_transform_edt(keep_lab_img)
ax2.imshow(dist_map, cmap='nipy_spectral')
ax2.set_title('Distance Map')
sns.heatmap(dist_map[::4, ::4],
annot=True,
fmt="1.0f",cmap='nipy_spectral',
ax=ax3,cbar=False,
vmin=0,vmax=dist_map.max(),
annot_kws={"size": 10})
ax3.set_title('Distance Map');
Skeletonization: Ridges¶
By using the Laplacian filter as an approximate for the derivative operator which finds the values which high local gradients.
$$ \nabla I_{d}(x,y) = (\frac{\delta^2}{\delta x^2}+\frac{\delta^2}{\delta y^2})I_d \approx \underbrace{\begin{bmatrix} -1 & -1 & -1 \\ -1 & 8 & -1 \\ -1 & -1 & -1 \end{bmatrix}}_{\textrm{Laplacian Kernel}} \otimes I_d(x,y) $$Creating the skeleton¶
We can locate the local maxima of the structure by setting a minimum surface distance
$$I_d(x,y)>MIN_{DIST}$$and combining it with a minimum slope value
$$\nabla I_d(x,y) > MIN_{SLOPE}$$Thresholds¶
Harking back to our earlier lectures, this can be seen as a threshold on a feature vector representation of the entire dataset.
- We first make the dataset into a tuple
In [65]:
# for finding the medial axis and making skeletons
from skimage.morphology import medial_axis
# for just the skeleton code
from skimage.morphology import skeletonize, skeletonize_3d
from skimage.filters import laplace
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 7))
ax1.imshow(dist_map, cmap='nipy_spectral');ax1.set_title('Distance Map')
ax2.imshow(laplace(dist_map), cmap='RdBu'); ax2.set_title('Laplacian of Distance')
# we use medial axis since it is cleaner
skel = medial_axis(keep_lab_img, return_distance=False) ; ax3.imshow(skel, cmap='gray'); ax3.set_title('Distance Map Ridge');
Skeletonization with Morphological thinning¶
From scikit-image documentation (http://scikit-image.org/docs/dev/auto_examples/edges/plot_skeleton.html)
Morphological thinning, implemented in the thin function, works on the same principle as skeletonize:
- remove pixels from the borders at each iteration until none can be removed without altering the connectivity.
- The different rules of removal can speed up skeletonization and result in different final skeletons.
The thin function also takes an optional max_iter keyword argument to limit the number of thinning iterations, and thus produce a relatively thicker skeleton.```
We can use this to thin the tiny junk elements first then erode, then perform the full skeletonization
Try morphological thinning¶
In [66]:
from skimage.morphology import thin, erosion
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(20, 7))
ax1.imshow(keep_lab_img, cmap="gray"); ax1.set_title('Segmentation')
thin_image = thin(keep_lab_img, max_iter=1)
ax2.imshow(thin_image, cmap="gray"); ax2.set_title('Morphologically Thinned')
er_thin_image = opening(thin_image, disk(1))
er_thin_image = label(er_thin_image) == 1
ax3.imshow(er_thin_image, cmap="gray"); ax3.set_title('Opened')
opened_skeleton = medial_axis(er_thin_image, return_distance=False)
ax4.imshow(opened_skeleton, cmap="gray"); ax4.set_title('Thinned/Opened Skeleton');
Still overgrown¶
The skeleton is still problematic for us and so we require some additional improvements to get a perfect skeleton
Skeleton: Junction Overview¶
With the skeleton which is ideally one voxel thick, we can characterize the junctions in the system by looking at the neighborhood of each voxel.
In [67]:
from scipy.ndimage import convolve
fig, ( ax1, ax3,ax2) = plt.subplots(1, 3, figsize=(18, 7))
ax1.imshow(opened_skeleton)
ax1.set_title('Skeleton')
neighbor_conv = convolve(opened_skeleton.astype(int), np.ones((3, 3)))
ax3.imshow(neighbor_conv,
cmap='nipy_spectral',
vmin=1, vmax=4,
interpolation='none' )
neighbor_conv[~opened_skeleton] = 0
j_img = ax2.imshow(neighbor_conv,
cmap='nipy_spectral',
vmin=1, vmax=4,
interpolation='none')
plt.colorbar(j_img)
ax2.set_title('Junction Count')
Out[67]:
Text(0.5, 1.0, 'Junction Count')
In [68]:
fig, ax1 = plt.subplots(1, 1, figsize=(8, 12))
n_crop = neighbor_conv[50:90, 40:60]
sns.heatmap(n_crop,
annot=True,
fmt="d",
cmap='nipy_spectral',
ax=ax1,
cbar=True,
vmin=0,
vmax=n_crop.max(),
annot_kws={"size": 10});
In [69]:
junc_types = np.unique(neighbor_conv[neighbor_conv > 0])
fig, m_axs = plt.subplots(1, len(junc_types), figsize=(20, 7))
for i, c_ax in zip(junc_types, m_axs):
c_ax.imshow(neighbor_conv == i, interpolation='none')
c_ax.set_title('Count == {}'.format(i))
Smarter kernel coding¶
Using a kernel like this: $$ j_4=\begin{array}{|c|c|c|} \hline & 1 & \\ \hline 2 & 256 & 4 \\ \hline & 8 & \\ \hline\end{array} \qquad j_8=\begin{array}{|c|c|c|} \hline 1 & 2 & 4 \\ \hline 8 & 256 & 16 \\ \hline 32 & 64 & 128 \\ \hline\end{array} $$
The neighborhood is uniquely coded.
Coding pixels with values $x$:
- $x = 0$ are background
- $x < 256$ touch the skeleton
- $x > 256$ are on the skeleton
The number of branches and their orientation can be encoded by counting bit flips or using a LUT.
In [70]:
fig, (ax1,ax2,ax3) = plt.subplots(1, 3, figsize=(15, 9))
n_crop = neighbor_conv[50:90, 40:60]
neighbor_j4 = convolve(opened_skeleton[50:90, 40:60].astype(int), np.array([[0,1,0],[2,16,4],[0,8,0]]))
neighbor_j8 = convolve(opened_skeleton[50:90, 40:60].astype(int), np.array([[1,2,4],[8,256,16],[32,64,128]]))
sns.heatmap(n_crop, annot=True, fmt="d",
cmap='nipy_spectral', ax=ax1, cbar=True,
vmin=0, vmax=n_crop.max(), annot_kws={"size": 6}), ax1.set_title('Plain convolution')
sns.heatmap(neighbor_j4, annot=True, fmt="d",
cmap='nipy_spectral', ax=ax2, cbar=True,
vmin=0, vmax=neighbor_j4.max(), annot_kws={"size": 6}), ax2.set_title('')
sns.heatmap(neighbor_j8, annot=True, fmt="d",
cmap='nipy_spectral', ax=ax3, cbar=True,
vmin=0, vmax=neighbor_j8.max(), annot_kws={"size": 6});
Dedicated pruning algorithms¶
- Ideally model-based
- Minimum branch length (using component labeling on the Count==3)
- Minimum branch width (using the distance map values)
In [71]:
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 6))
lab_seg = label(neighbor_conv == 3)
ax1.imshow(lab_seg, cmap='gist_earth')
ax1.set_title('Segment ID')
ax2.hist(lab_seg[lab_seg > 0])
ax2.set_title('Segment Length')
label_length_img = np.zeros_like(lab_seg)
for i in np.unique(lab_seg[lab_seg > 0]):
label_length_img[lab_seg == i] = np.sum(lab_seg == i)
ll_ax = ax3.imshow(label_length_img, cmap='jet')
ax3.set_title('Segment Length'); plt.colorbar(ll_ax);
In [72]:
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 7))
length_skeleton = (label_length_img > 5) + \
(neighbor_conv == 2) + \
(neighbor_conv > 3)
ax1.imshow(im_crop)
ax2.imshow(length_skeleton)
ax3.imshow(length_skeleton)
ax3.scatter(mk_crop['x'], mk_crop['y'], s=mk_crop['width'],
alpha=0.25, label='Ground Truth')
ax3.legend();
Analyzing maximum segment width¶
In [73]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 6))
label_width_img = np.zeros_like(lab_seg)
for i in np.unique(lab_seg[lab_seg > 0]):
label_width_img[lab_seg == i] = np.max(dist_map[lab_seg == i])
ax1.hist(label_width_img[label_width_img > 0])
ax1.set_title('Segment Maximum Width')
ll_ax = ax2.imshow(label_width_img, cmap='jet')
ax2.set_title('Segment Maximum Width'); plt.colorbar(ll_ax);
Pruning using structure width¶
In [92]:
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(20, 7))
width_skeleton = (label_width_img > 4.5) \
+ (neighbor_conv == 2) \
+ (neighbor_conv > 3)
width_skeleton = label(width_skeleton) == 1
ax1.imshow(im_crop)
ax2.imshow(width_skeleton), ax2.set_title('Pruned skeleton')
ax3.imshow(width_skeleton)
ax3.scatter(mk_crop['x'], mk_crop['y'], s=mk_crop['width'],
alpha=0.25, label='Ground Truth')
ax3.legend();
ax4.imshow(medial_axis(keep_lab_img, return_distance=False)), ax4.set_title('Medial axis skeleton')
Out[92]:
(<matplotlib.image.AxesImage at 0x1eba0c20b48>, Text(0.5, 1.0, 'Medial axis skeleton'))
Establish Topology¶
From the cleaned, pruned skeleton we can start to establish topology.
Using the same criteria as before we can break down the image into
- segments,
- junctions,
- and end-points
In [75]:
ws_neighbors = convolve(width_skeleton.astype(
int), np.ones((3, 3)), mode='constant', cval=0)
ws_neighbors[~width_skeleton] = 0
fig, (ax1) = plt.subplots(1, 1, figsize=(5,10), dpi=100)
ax1.imshow(im_crop)
j_name = {1: 'dangling point', 2: 'end-point',
3: 'segment', 4: 'junction', 5: 'super-junction'}
for j_count in np.unique(ws_neighbors[ws_neighbors > 0]):
y_c, x_c = np.where(ws_neighbors == j_count)
ax1.plot(x_c, y_c, 's',
label=j_name.get(j_count, 'unknown'),
markersize=5)
leg = ax1.legend(shadow=True, fancybox=True, frameon=True)
Getting Topology in Image Space¶
We want to determine which nodes are directly connected in this image so we can extract a graph. If we take a simple case of two nodes connected by one edge and the bottom node connected to another edge going nowhere.
$$ \begin{bmatrix} n & 0 & 0 & 0 \\ 0 & e & 0 & 0 \\ 0 & 0 & n & e \end{bmatrix} $$We can use component labeling to identify each node and each edge uniquely
Node Labels¶
$$ N_{lab} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 2 & 0 \end{bmatrix} $$Edge Labels¶
$$E_{lab} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 2 \end{bmatrix} $$We can then use a dilation operation on the nodes and the edges to see which overlap
In [78]:
from skimage.morphology import dilation
n_img = np.zeros((3, 4))
e_img = np.zeros_like(n_img)
n_img[0, 0] = 1
e_img[1, 1] = 1
n_img[2, 2] = 1
e_img[2, 3] = 1
fig, ((ax1, ax3, ax5), (ax2, ax4, ax6)) = plt.subplots(2, 3, figsize=(20,8))
ax1.imshow(n_img)
ax1.set_title('Nodes')
ax2.imshow(e_img)
ax2.set_title('Edges')
# labeling
n_labs = label(n_img)
sns.heatmap(n_labs, annot=True, fmt="d", ax=ax3, cbar=False); ax3.set_title('Node Labels')
e_labs = label(e_img)
sns.heatmap(e_labs, annot=True, fmt="d", ax=ax4, cbar=False); ax4.set_title('Edge Labels')
# growing
n_grow_1 = dilation(n_labs == 1, np.ones((3, 3)))
sns.heatmap(n_grow_1, annot=True, fmt="d", ax=ax5, cbar=False);
ax5.set_title('Grow First\n{} {}'.format('Edges Found:', [x for x in np.unique(e_labs[n_grow_1 > 0]) if x > 0]))
n_grow_2 = dilation(n_labs == 2, np.ones((3, 3)))
sns.heatmap(n_grow_2, annot=True, fmt="d", ax=ax6, cbar=False)
ax6.set_title('Grow Second\n{} {}'.format('Edges Found:', [x for x in np.unique(e_labs[n_grow_2 > 0]) if x > 0]));
In [77]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))
node_id_image = label((ws_neighbors > 3) | (ws_neighbors == 2))
edge_id_image = label(ws_neighbors == 3)
ax1.imshow(im_crop)
node_dict = {}
for c_node in np.unique(node_id_image[node_id_image > 0]):
y_n, x_n = np.where(node_id_image == c_node)
node_dict[c_node] = {'x': np.mean(x_n),
'y': np.mean(y_n),
'width': np.mean(dist_map[node_id_image == c_node])}
ax1.plot(np.mean(x_n), np.mean(y_n), 'rs')
edge_dict = {}
edge_matrix = np.eye(len(node_dict)+1)
for c_edge in np.unique(edge_id_image[edge_id_image > 0]):
edge_grow_mask = dilation(edge_id_image == c_edge, np.ones((3, 3)))
v_nodes = np.unique(node_id_image[edge_grow_mask > 0])
v_nodes = [v for v in v_nodes if v > 0]
print('Edge', c_edge, 'connects', v_nodes)
if len(v_nodes) == 2:
edge_dict[c_edge] = {'start': v_nodes[0],
'end': v_nodes[-1],
'length': np.sum(edge_id_image == c_edge),
'euclidean_distance': np.sqrt(np.square(node_dict[v_nodes[0]]['x'] -
node_dict[v_nodes[-1]]['x']) +
np.square(node_dict[v_nodes[0]]['y'] -
node_dict[v_nodes[-1]]['y'])
),
'max_width': np.max(dist_map[edge_id_image == c_edge]),
'mean_width': np.mean(dist_map[edge_id_image == c_edge])}
edge_matrix[v_nodes[0], v_nodes[-1]] = np.sum(edge_id_image == c_edge)
edge_matrix[v_nodes[-1], v_nodes[0]] = np.sum(edge_id_image == c_edge)
s_node = node_dict[v_nodes[0]]
e_node = node_dict[v_nodes[-1]]
ax1.plot([s_node['x'], e_node['x']],
[s_node['y'], e_node['y']], 'b-', linewidth=np.mean(dist_map[edge_id_image == c_edge]), alpha=0.5)
ax2.matshow(edge_matrix, cmap='viridis')
ax2.set_title('Connectivity Matrix'); ax2.set_xlabel('Node A'); ax2.set_ylabel('Node B');
Edge 1 connects [1, 3] Edge 2 connects [2, 3] Edge 3 connects [3, 4] Edge 4 connects [3, 6] Edge 5 connects [5, 6] Edge 6 connects [6, 7] Edge 7 connects [6, 8] Edge 8 connects [8, 9] Edge 9 connects [8, 10] Edge 10 connects [9] Edge 11 connects [10, 11] Edge 12 connects [10, 12] Edge 13 connects [11] Edge 14 connects [11, 14] Edge 15 connects [11, 13] Edge 16 connects [14] Edge 17 connects [14, 15] Edge 18 connects [14, 16] Edge 19 connects [16, 17] Edge 20 connects [16, 17] Edge 21 connects [17, 18] Edge 22 connects [18] Edge 23 connects [18, 19] Edge 24 connects [19, 20] Edge 25 connects [19, 21] Edge 26 connects [21] Edge 27 connects [21, 23] Edge 28 connects [22, 23] Edge 29 connects [23, 24] Edge 30 connects [24] Edge 31 connects [24] Edge 32 connects [24, 25] Edge 33 connects [25, 26] Edge 34 connects [26] Edge 35 connects [26, 27] Edge 36 connects [25, 28]
Skeleton: Tortuosity¶
One of the more interesting ones in material science is called tortuosity and it is defined as the ratio between the arc-length of a segment and the distance between its starting and ending points. $$ \tau = \frac{L}{C} $$
A high degree of tortuosity indicates that the network is convoluted and is important when estimating or predicting flow rates. Specifically
- in geology it is an indication that diffusion and fluid transport will occur more slowly
- in analytical chemistry it is utilized to perform size exclusion chromatography
- in vascular tissue it can be a sign of pathology.
In [79]:
fig, (ax1) = plt.subplots(1, 1, figsize=(8, 8))
ax1.imshow(im_crop)
for _, d_values in edge_dict.items():
v_nodes = [d_values['start'], d_values['end']]
print('Tortuousity: %2.2f' %
(d_values['length']/d_values['euclidean_distance']))
s_node = node_dict[v_nodes[0]]
e_node = node_dict[v_nodes[-1]]
ax1.plot([s_node['x'], e_node['x']],
[s_node['y'], e_node['y']], 'b-',
linewidth=5, alpha=d_values['length']/d_values['euclidean_distance'])
Tortuousity: 0.75 Tortuousity: 0.66 Tortuousity: 0.80 Tortuousity: 0.76 Tortuousity: 0.92 Tortuousity: 0.78 Tortuousity: 0.87 Tortuousity: 0.95 Tortuousity: 0.88 Tortuousity: 0.61 Tortuousity: 0.71 Tortuousity: 0.23 Tortuousity: 0.67 Tortuousity: 0.67 Tortuousity: 0.83 Tortuousity: 0.50 Tortuousity: 0.50 Tortuousity: 0.77 Tortuousity: 0.89 Tortuousity: 0.86 Tortuousity: 0.31 Tortuousity: 0.84 Tortuousity: 0.54 Tortuousity: 0.22 Tortuousity: 0.67 Tortuousity: 0.89 Tortuousity: 0.91 Tortuousity: 0.94
In [80]:
import networkx as nx
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))
G = nx.Graph()
for k, v in node_dict.items():
G.add_node(k, weight=v['width'])
for k, v in edge_dict.items():
G.add_edge(v['start'], v['end'], **v)
nx.draw_spring(G, ax=ax1, with_labels=True,
node_color=[node_dict[k]['width']
for k in sorted(node_dict.keys())],
node_size=800,
cmap=plt.cm.autumn,
edge_color=[G.edges[k]['length'] for k in list(G.edges.keys())],
width=[2*G.edges[k]['max_width'] for k in list(G.edges.keys())],
edge_cmap=plt.cm.Greens)
ax1.set_title('Randomly Organized Graph')
ax2.imshow(im_crop)
nx.draw(G,
pos={k: (v['x'], v['y']) for k, v in node_dict.items()},
ax=ax2,
node_color=[node_dict[k]['width'] for k in sorted(node_dict.keys())],
node_size=50,
cmap=plt.cm.autumn,
edge_color=[G.edges[k]['length'] for k in list(G.edges.keys())],
width=[2*G.edges[k]['max_width'] for k in list(G.edges.keys())],
edge_cmap=plt.cm.Blues,
alpha=0.5,
with_labels=False)
Graph Analysis¶
Once the data has been represented in a graph form, we can begin to analyze some of graph aspects of it, like the degree and connectivity plots.
In [94]:
degree_sequence = sorted([d for n, d in G.degree()],
reverse=True) # degree sequence
plt.hist(degree_sequence, bins=np.arange(10)), plt.xlabel('Number of connected edges'), plt.ylabel('Number of nodes');
Skeletons going 3D¶
- Object topology must still be preserved
- More complex neighborhoods to analyze
Load volume¶
The data is zipped in the repos, unzip before continuing
In [83]:
root = np.load('../common/data/Cropped_prediction_8bit.npy')
fig,(ax1,ax2,ax3) = plt.subplots(1,3,figsize=(20,5))
ax1.imshow(root[:,:,700]); ax1.set_title('XY')
ax2.imshow(root[:,220,:]); ax2.set_title('XZ')
ax3.imshow(root[220,:,:]); ax3.set_title('YZ');
Let look at the projection instead¶
In [88]:
plt.figure(figsize=(15,4))
plt.imshow(root.mean(axis=0)), plt.title('Projection of the segmented image');
Create the 3D skeleton¶
In [85]:
skel = skeletonize_3d(root)
In [89]:
plt.figure(figsize=(15,4))
plt.imshow(skel.max(axis=0), cmap="bone"), plt.title('Projection of the skeleton');
Segmenting touching items¶
Watershed¶
Watershed is a method for segmenting objects without using component labeling.
- It utilizes the shape of structures to find objects
- From the distance map we can make out substructures with our eyes
- But how to we find them?!
We use a sample image now from the Datascience Bowl 2018 from Kaggle. The challenge is to identify nuclei in histology images to eventually find cancer better. The winner tweeted about the solution here
In [100]:
from skimage.filters import threshold_otsu
from skimage.color import rgb2hsv
import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread
%matplotlib inline
rgb_img = imread("../common/figures/dsb_sample/slide.png")[:, :, :3]
gt_labs = imread("../common/figures/dsb_sample/labels.png")
bw_img = rgb2hsv(rgb_img)[:, :, 2]
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(20, 6), dpi=100)
ax1.imshow(rgb_img, cmap='bone')
ax2.imshow(bw_img, cmap='bone'),ax2.set_title('Gray Scale')
ax3.imshow(bw_img < threshold_otsu(bw_img), cmap='bone'), ax3.set_title('Segmentation')
ax4.imshow(gt_labs, cmap='tab20c'),ax4.set_title('Labeled cells');
Lets try component labeling¶
In [101]:
from skimage.morphology import label
import seaborn as sns
fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(20, 6), dpi=100)
bw_roi = bw_img[75:110:2, 125:150:2]
ax1.imshow(bw_roi, cmap='bone'); ax1.set_title('Gray Scale')
bw_roi_seg = bw_roi < threshold_otsu(bw_img)
sns.heatmap(bw_roi_seg, annot=True, fmt="d",
ax=ax2, cbar=False, cmap='gist_earth'); ax2.set_title('Segmentation');
bw_roi_label = label(bw_roi_seg)
sns.heatmap(bw_roi_label, annot=True, fmt="d",
ax=ax3, cbar=False, cmap='gist_earth'); ax3.set_title('Labels')
sns.heatmap(gt_labs[75:110:2, 125:150:2], annot=True,
fmt="d", ax=ax4, cbar=False, cmap='gist_earth'); ax4.set_title('Ground Truth');
Watershed: Flowing Downhill¶
We can imagine watershed as waterflowing down hill into basins. The topology in this case is given by the distance map
In [102]:
from scipy.ndimage import distance_transform_edt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(15, 10))
ax = fig.gca(projection='3d')
bw_roi_dmap = distance_transform_edt(bw_roi_seg) # The distance map needed for the segmentation
# Plot the surface.
t_xx, t_yy = np.meshgrid(np.arange(bw_roi_dmap.shape[1]),np.arange(bw_roi_dmap.shape[0]))
surf = ax.plot_surface(t_xx, t_yy,-1*bw_roi_dmap, cmap="terrain",linewidth=0.25, antialiased=True)
# Customize the z axis.
ax.view_init(60, 20)
# Add a color bar which maps values to colors.
fig.colorbar(surf, shrink=0.5); ax.set_title('Elevation map',fontsize=16);
Preparations for watershed segmentation¶
We need to create a seed image
- Segment image
- Compute distance map
- Identify local maxima
In [103]:
from skimage.feature import peak_local_max
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 6), dpi=100)
sns.heatmap(bw_roi_seg, annot=True, fmt="d",
ax=ax1, cbar=False, cmap='gist_earth'); ax1.set_title('Segmentation')
sns.heatmap(bw_roi_dmap, annot=True, fmt="1.0f",ax=ax2, cbar=False, cmap='viridis'),ax2.set_title('Distance Map');
roi_local_maxi = peak_local_max(bw_roi_dmap, indices=False, footprint=np.ones((3, 3)),
labels=bw_roi_seg, exclude_border=False)
labeled_maxi = label(roi_local_maxi)
sns.heatmap(labeled_maxi, annot=True, fmt="1.0f", ax=ax3, cbar=False, cmap='gist_earth'),ax3.set_title('Local Maxima');
Run watershed¶
In [104]:
from skimage.morphology import watershed
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 6), dpi=100)
sns.heatmap(labeled_maxi, annot=True, fmt="1.0f",
ax=ax1, cbar=False, cmap='gist_earth'); ax1.set_title('Local Maxima')
ws_labels = watershed(-bw_roi_dmap, labeled_maxi, mask=bw_roi_seg)
sns.heatmap(ws_labels, annot=True, fmt="d",
ax=ax2, cbar=False, cmap='gist_earth'); ax2.set_title('Watershed')
sns.heatmap(gt_labs[75:110:2, 125:150:2], annot=True,
fmt="d", ax=ax3, cbar=False, cmap='gist_earth'); ax3.set_title('Ground Truth');
Removing too small elements¶
One of the components (Label=2) is too small!
We can remove it by deleting unwanted seeds.
- seed belonging to the bottom 10 percentile of areas
- rerunning watershed
In [105]:
label_area_dict = {i: np.sum(ws_labels == i)
for i in np.unique(ws_labels[ws_labels > 0])}
clean_label_maxi = labeled_maxi.copy()
area_cutoff = np.percentile(list(label_area_dict.values()), 10)
print('!0% cutoff', area_cutoff)
for i, k in label_area_dict.items():
print('Label: ', i, 'Area:', k, 'Keep:', k > area_cutoff)
if k <= area_cutoff:
clean_label_maxi[clean_label_maxi == i] = 0
!0% cutoff 11.2 Label: 1 Area: 82 Keep: True Label: 2 Area: 7 Keep: False Label: 3 Area: 59 Keep: True Label: 4 Area: 21 Keep: True
In [106]:
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 6), dpi=100)
sns.heatmap(clean_label_maxi, annot=True, fmt="1.0f",
ax=ax1, cbar=False, cmap='gist_earth'); ax1.set_title('Local Maxima')
ws_labels = watershed(-bw_roi_dmap, clean_label_maxi, mask=bw_roi_seg)
sns.heatmap(ws_labels, annot=True, fmt="d",
ax=ax2, cbar=False, cmap='gist_earth'); ax2.set_title('Watershed')
sns.heatmap(gt_labs[75:110:2, 125:150:2], annot=True,
fmt="d", ax=ax3, cbar=False, cmap='gist_earth');ax3.set_title('Ground Truth');
Scaling back up¶
Now we can perform the operation on the whole image and see how the results look
In [107]:
from skimage.morphology import opening, disk
from skimage.segmentation import mark_boundaries
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 6), dpi=150)
ax1.imshow(rgb_img, cmap='bone')
bw_seg_img = opening(bw_img < threshold_otsu(bw_img), disk(3))
bw_dmap = distance_transform_edt(bw_seg_img)
bw_peaks = label(peak_local_max(bw_dmap, indices=False, footprint=np.ones((3, 3)),
labels=bw_seg_img, exclude_border=True))
ws_labels = watershed(-bw_dmap, bw_peaks, mask=bw_seg_img)
label_area_dict = {i: np.sum(ws_labels == i)
for i in np.unique(ws_labels[ws_labels > 0])}
clean_label_maxi = bw_peaks.copy()
lab_areas = list(label_area_dict.values())
area_cutoff = np.percentile(lab_areas, 20)
print('10% cutoff', area_cutoff, 'Removed', np.sum(
np.array(lab_areas) < area_cutoff), 'components')
for i, k in label_area_dict.items():
if k <= area_cutoff:
clean_label_maxi[clean_label_maxi == i] = 0
ws_labels = watershed(-bw_dmap, clean_label_maxi, mask=bw_seg_img)
ax2.imshow(mark_boundaries(label_img=ws_labels,
image=rgb_img, color=(0, 1, 0))); ax2.set_title('Watershed')
ax3.imshow(mark_boundaries(label_img=gt_labs, image=rgb_img, color=(0, 1, 0))); ax3.set_title('Ground Truth');
plt.tight_layout(); fig.savefig('ws_full.png',dpi=300);
10% cutoff 34.8 Removed 24 components
Why bad performance?¶
- Too simple thresholding method was used (Otsu)
- Irregular shapes produce too many local maxima in the elevation map
Solution¶
- Find better segmentation approach than only histogram
- Use morphological algorithms like max-impose Soille, 1999
Battery Example¶
We use an example from the Laboratory for Nanoelectronics at ETH Zurich. The datasets are x-ray tomography images of the battery micro- and nanostructures. As the papers below document, a substantial amount of image processing is required to extract meaningful physical and chemical values from these images.
Acknowledgements¶
The relevant publications which also contain links to the full collection of datasets.
Goal¶
The goal is to
- Segment and quantify the relevant structures
- Find out what changes occur between 0 and 2000 bar of pressure.
In [108]:
from scipy.ndimage import binary_fill_holes
from skimage.morphology import opening, closing, disk
from skimage.filters import threshold_otsu
import numpy as np
import matplotlib.pyplot as plt
from skimage.io import imread
%matplotlib inline
bw_img = imread("../common/data/NMC_90wt_2000bar_115.tif")[:, :, 0]
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 10), dpi=120)
ax1.imshow(bw_img, cmap='bone'), ax1.set_title('Gray Scale')
thresh_img = bw_img > threshold_otsu(bw_img)
ax2.imshow(thresh_img, cmap='bone'), ax2.set_title('Segmentation (Otsu)')
bw_seg_img = closing(
closing(
opening(thresh_img, disk(3)),
disk(1)
), disk(1)
)
ax3.imshow(bw_seg_img, cmap='bone'); ax3.set_title('Clean Segments after closing');
Let the water flow...¶
In [109]:
from scipy.ndimage import distance_transform_edt
from skimage.morphology import label
from skimage.feature import peak_local_max
from skimage.segmentation import mark_boundaries
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 8), dpi=100)
bw_dmap = distance_transform_edt(bw_seg_img)
ax1.imshow(bw_dmap, cmap='nipy_spectral'), ax1.set_title('Distance map')
bw_peaks = label(peak_local_max(bw_dmap, indices=False, footprint=np.ones((3, 3)),
labels=bw_seg_img, exclude_border=True))
ws_labels = watershed(-bw_dmap, bw_peaks, mask=bw_seg_img)
ax2.imshow(ws_labels, cmap='gist_earth'), ax2.set_title('Watershed labels')
# find boundaries
ax3.imshow(mark_boundaries(label_img=ws_labels, image=bw_img)); ax3.set_title('Watershed boundaries');
Too much detail...¶
In [110]:
def get_roi(x): return x[0:200, 200:450]
im_crop = get_roi(bw_img)
dist_map = get_roi(bw_dmap)
node_id_image = get_roi(ws_labels)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 8))
ax1.imshow(im_crop), ax1.set_title('Cropped image')
ax2.imshow(dist_map), ax2.set_title('Distance map')
#ax3.imshow(node_id_image,cmap='gist_earth'), ax3.set_title('Watershed labels');
ax3.imshow(node_id_image,cmap='flag'), ax3.set_title('Watershed labels');
Representing as a Graph¶
Here we can change the representation from a number of random labels to a graph
In [111]:
from skimage.morphology import dilation
from skimage.measure import perimeter
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))
ax1.imshow(im_crop)
node_dict = {}
for c_node in np.unique(node_id_image[node_id_image > 0]):
y_n, x_n = np.where(node_id_image == c_node)
node_dict[c_node] = {'x': np.mean(x_n),
'y': np.mean(y_n),
'width': np.mean(dist_map[node_id_image == c_node]),
'perimeter': perimeter(node_id_image == c_node)}
ax1.plot(np.mean(x_n), np.mean(y_n), 'rs')
edge_dict = {}
for i in node_dict.keys():
i_grow = dilation(node_id_image == i, np.ones((3, 3)))
for j in node_dict.keys():
if i < j:
j_grow = dilation(node_id_image == j, np.ones((3, 3)))
interface_length = np.sum(i_grow & j_grow)
if interface_length > 0:
v_nodes = [i, j]
edge_dict[(i, j)] = {'start': v_nodes[0],
'start_perimeter': node_dict[v_nodes[0]]['perimeter'],
'end_perimeter': node_dict[v_nodes[-1]]['perimeter'],
'end': v_nodes[-1],
'interface_length': interface_length,
'euclidean_distance': np.sqrt(np.square(node_dict[v_nodes[0]]['x'] -
node_dict[v_nodes[-1]]['x']) +
np.square(node_dict[v_nodes[0]]['y'] -
node_dict[v_nodes[-1]]['y'])
),
'max_width': np.max(dist_map[i_grow & j_grow]),
'mean_width': np.mean(dist_map[i_grow & j_grow])}
s_node = node_dict[v_nodes[0]]
e_node = node_dict[v_nodes[-1]]
ax1.plot([s_node['x'], e_node['x']],
[s_node['y'], e_node['y']], 'b-',
linewidth=np.max(dist_map[i_grow & j_grow]), alpha=0.5)
ax2.imshow(mark_boundaries(label_img=node_id_image, image=im_crop))
ax2.set_title('Borders');
In [112]:
import pandas as pd
edge_df = pd.DataFrame(list(edge_dict.values()))
edge_df.head(5)
Out[112]:
| start | start_perimeter | end_perimeter | end | interface_length | euclidean_distance | max_width | mean_width | |
|---|---|---|---|---|---|---|---|---|
| 0 | 6 | 31.656854 | 25.621320 | 10 | 32 | 3.517456 | 6.403124 | 3.191413 |
| 1 | 6 | 31.656854 | 27.656854 | 18 | 1 | 6.265257 | 6.324555 | 6.324555 |
| 2 | 10 | 25.621320 | 27.656854 | 18 | 30 | 2.748111 | 6.324555 | 3.217437 |
| 3 | 18 | 27.656854 | 19.621320 | 27 | 28 | 2.612361 | 5.656854 | 2.666442 |
| 4 | 18 | 27.656854 | 25.313708 | 33 | 1 | 6.368174 | 5.000000 | 5.000000 |
Combine split electrodes¶
Here we combine split electrodes by using a cutoff on the ratio of the interface length to the start and end perimeters
In [113]:
delete_edges = edge_df.query(
'interface_length>0.33*(start_perimeter+end_perimeter)')
print('Found', delete_edges.shape[0], '/', edge_df.shape[0], 'edges to delete')
delete_edges.head(5)
Found 76 / 203 edges to delete
Out[113]:
| start | start_perimeter | end_perimeter | end | interface_length | euclidean_distance | max_width | mean_width | |
|---|---|---|---|---|---|---|---|---|
| 0 | 6 | 31.656854 | 25.621320 | 10 | 32 | 3.517456 | 6.403124 | 3.191413 |
| 2 | 10 | 25.621320 | 27.656854 | 18 | 30 | 2.748111 | 6.324555 | 3.217437 |
| 3 | 18 | 27.656854 | 19.621320 | 27 | 28 | 2.612361 | 5.656854 | 2.666442 |
| 6 | 27 | 19.621320 | 25.313708 | 33 | 24 | 3.874383 | 5.000000 | 2.409123 |
| 7 | 28 | 12.414214 | 28.727922 | 29 | 16 | 5.028904 | 3.000000 | 1.404919 |
In [117]:
node_id_image = get_roi(ws_labels)
for _ in range(3):
# since some mappings might be multistep
for _, c_row in delete_edges.iterrows():
node_id_image[node_id_image == c_row['end']] = c_row['start']
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))
ax1.imshow(im_crop)
node_dict = {}
for c_node in np.unique(node_id_image[node_id_image > 0]):
y_n, x_n = np.where(node_id_image == c_node)
node_dict[c_node] = {'x': np.mean(x_n),
'y': np.mean(y_n),
'width': np.mean(dist_map[node_id_image == c_node]),
'perimeter': perimeter(node_id_image == c_node)}
ax1.plot(np.mean(x_n), np.mean(y_n), 'rs')
edge_dict = {}
for i in node_dict.keys():
i_grow = dilation(node_id_image == i, np.ones((3, 3)))
for j in node_dict.keys():
if i < j:
j_grow = dilation(node_id_image == j, np.ones((3, 3)))
interface_length = np.sum(i_grow & j_grow)
if interface_length > 0:
v_nodes = [i, j]
edge_dict[(i, j)] = {'start': v_nodes[0],
'start_perimeter': node_dict[v_nodes[0]]['perimeter'],
'end_perimeter': node_dict[v_nodes[-1]]['perimeter'],
'end': v_nodes[-1],
'interface_length': interface_length,
'euclidean_distance': np.sqrt(np.square(node_dict[v_nodes[0]]['x'] -
node_dict[v_nodes[-1]]['x']) +
np.square(node_dict[v_nodes[0]]['y'] -
node_dict[v_nodes[-1]]['y'])
),
'max_width': np.max(dist_map[i_grow & j_grow]),
'mean_width': np.mean(dist_map[i_grow & j_grow])}
s_node = node_dict[v_nodes[0]]
e_node = node_dict[v_nodes[-1]]
ax1.plot([s_node['x'], e_node['x']],
[s_node['y'], e_node['y']], 'b-',
linewidth=np.max(dist_map[i_grow & j_grow]), alpha=0.5)
ax2.imshow(mark_boundaries(label_img=node_id_image, image=im_crop))
ax2.set_title('Borders');
In [115]:
import networkx as nx
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))
G = nx.Graph()
for k, v in node_dict.items():
G.add_node(k, weight=v['width'])
for k, v in edge_dict.items():
G.add_edge(v['start'], v['end'], **v)
nx.draw_shell(G, ax=ax1, with_labels=True,
node_color=[node_dict[k]['width']
for k in sorted(node_dict.keys())],
node_size=400,
cmap=plt.cm.autumn,
edge_color=[G.edges[k]['interface_length']
for k in list(G.edges.keys())],
width=[2*G.edges[k]['max_width'] for k in list(G.edges.keys())],
edge_cmap=plt.cm.Greens)
ax1.set_title('Randomly Organized Graph')
ax2.imshow(im_crop, cmap="viridis")
nx.draw(G,
pos={k: (v['x'], v['y']) for k, v in node_dict.items()},
ax=ax2,
node_color=[node_dict[k]['width'] for k in sorted(node_dict.keys())],
node_size=50,
cmap=plt.cm.Greens,
edge_color=[G.edges[k]['interface_length']
for k in list(G.edges.keys())],
width=[2*G.edges[k]['max_width'] for k in list(G.edges.keys())],
edge_cmap=plt.cm.autumn,
alpha=0.75,
with_labels=False)
In [116]:
degree_sequence = sorted([d for n, d in G.degree()],
reverse=True) # degree sequence
plt.hist(degree_sequence, bins=np.arange(10)), plt.xlabel('Number of edges connect to node'), plt.ylabel('Number of nodes');
In [ ]: