With inconsistent or every changing illumination it may not be possible to apply the same threshold to every image.
%matplotlib inline
from skimage.io import imread
import matplotlib.pyplot as plt
import numpy as np
cell_img = imread("../common/figures/Cell_Colony.jpg")/255.0
np.random.seed(2018)
m_cell_imgs = [cell_img+k for k in np.random.uniform(-0.25, 0.25, size = 8)]
fig, m_axs = plt.subplots(3, 3, figsize = (12, 12), dpi = 72)
ax1 = m_axs[0,0]
for i, (c_ax, c_img) in enumerate(zip(m_axs.flatten()[1:], m_cell_imgs)):
ax1.hist(c_img.ravel(), np.linspace(0, 1, 20),
label = '{}'.format(i), alpha = 0.5, histtype = 'step')
c_ax.imshow(c_img, cmap = 'bone', vmin = 0, vmax = 1)
c_ax.axis('off')
ax1.set_yscale('log', nonposy = 'clip')
ax1.set_ylabel('Pixel Count (Log)')
ax1.set_xlabel('Intensity')
ax1.legend();
A constant threshold of 0.6 for the different illuminations
fig, m_axs = plt.subplots(len(m_cell_imgs), 2,
figsize = (4*2, 4*len(m_cell_imgs)),
dpi = 72)
THRESH_VAL = 0.6
for i, ((ax1, ax2), c_img) in enumerate(zip(m_axs,
m_cell_imgs)):
ax1.hist(c_img.ravel()[c_img.ravel()>THRESH_VAL],
np.linspace(0, 1, 100),
label = 'Above Threshold',
alpha = 0.5)
ax1.hist(c_img.ravel()[c_img.ravel()<THRESH_VAL],
np.linspace(0, 1, 100),
label = 'Below Threshold',
alpha = 0.5)
ax1.set_yscale('log', nonposy = 'clip')
ax1.set_ylabel('Pixel Count (Log)')
ax1.set_xlabel('Intensity')
ax1.legend()
ax2.imshow(c_img<THRESH_VAL,
cmap = 'bone',
vmin = 0, vmax = 1)
ax2.set_title('Volume Fraction (%2.2f%%)' % (100*np.mean(c_img.ravel()<THRESH_VAL)))
ax2.axis('off')
is easy using a threshold and size criteria (we know how big the cells should be)
is much more difficult because the small channels having radii on the same order of the pixel size are obscured by partial volume effects and noise.
%matplotlib inline
from skimage.io import imread
import matplotlib.pyplot as plt
import numpy as np
cortex_img = imread("ext-figures/cortex.png")
np.random.seed(2018)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3,
figsize = (12, 4), dpi = 72)
ax1.imshow(cortex_img, cmap = 'bone')
thresh_vals = np.linspace(cortex_img.min(), cortex_img.max(), 4+2)[:-1]
out_img = np.zeros_like(cortex_img)
for i, (t_start, t_end) in enumerate(zip(thresh_vals, thresh_vals[1:])):
thresh_reg = (cortex_img>t_start) & (cortex_img<t_end)
ax2.hist(cortex_img.ravel()[thresh_reg.ravel()])
out_img[thresh_reg] = i
ax3.imshow(out_img, cmap = 'gist_earth');
Given that applying a threshold is such a common and signficant step, there have been many tools developed to automatically (unsupervised) perform it. A particularly important step in setups where images are rarely consistent such as outdoor imaging which has varying lighting (sun, clouds). The methods are based on several basic principles.
Just like we visually inspect a histogram an algorithm can examine the histogram and find local minimums between two peaks, maximum / minimum entropy and other factors
These look at the statistics of the thresheld image themselves (like entropy) to estimate the threshold
These search for a threshold which delivers the desired results in the final objects. For example if you know you have an image of cells and each cell is between 200-10000 pixels the algorithm runs thresholds until the objects are of the desired size
Taking a typical image of a bone slice, we can examine the variations in calcification density in the image
We can see in the histogram that there are two peaks, one at 0 (no absorption / air) and one at 0.5 (stronger absorption / bone)
%matplotlib inline
from skimage.io import imread
import matplotlib.pyplot as plt
import numpy as np
bone_img = imread("ext-figures/bonegfiltslice.png")/255.0
np.random.seed(2018)
fig, (ax1, ax2) = plt.subplots(1, 2,figsize = (12, 4), dpi = 72)
ax1.imshow(bone_img, cmap = 'bone')
ax1.axis('off')
ax2.hist(bone_img.ravel(), 50);
Search and minimize intra-class (within) variance $$\sigma^2_w(t)=\omega_{bg}(t)\sigma^2_{bg}(t)+\omega_{fg}(t)\sigma^2_{fg}(t)$$
There are many methods and they can be complicated to implement yourself. Both Fiji and scikit-image offers many of them as built in functions so you can automatically try all of them on your image
%matplotlib inline
from skimage.io import imread
import matplotlib.pyplot as plt
import numpy as np
bone_img = imread("ext-figures/bonegfiltslice.png")/255.0
from skimage.filters import try_all_threshold
fig, ax = try_all_threshold(bone_img, figsize=(10, 8), verbose=True)
skimage.filters.thresholding.threshold_isodata skimage.filters.thresholding.threshold_li skimage.filters.thresholding.threshold_mean skimage.filters.thresholding.threshold_minimum skimage.filters.thresholding.threshold_otsu skimage.filters.thresholding.threshold_triangle skimage.filters.thresholding.threshold_yen
While an incredibly useful tool, there are many potential pitfalls to these automated techniques.
These methods are very sensitive to the distribution of pixels in your image and may work really well on images with equal amounts of each phase but work horribly on images which have very high amounts of one phase compared to the others
These methods are sensitive to noise and a large noise content in the image can change statistics like entropy significantly.
These methods are inherently biased by the expectations you have. If you want to find objects between 200 and 1000 pixels you will, they just might not be anything meaningful.
Imaging science rarely represents the ideal world and will never be 100% perfect. At some point we need to write our master's thesis, defend, or publish a paper. These are approaches for more qualitative assessment we will later cover how to do this a bit more robustly with quantitative approaches
One approach is to try and simulate everything (including noise) as well as possible and to apply these techniques to many realizations of the same image and qualitatively keep track of how many of the results accurately identify your phase or not. Hint: >95% seems to convince most biologists
Apply the methods to each sample and keep track of which threshold was used for each one. Go back and apply each threshold to each sample in the image and keep track of how many of them are correct enough to be used for further study.
Come up with the worst-case scenario (noise, misalignment, etc) and assess how inacceptable the results are. Then try to estimate the quartiles range (75% - 25% of images).
For some images a single threshold does not work
%matplotlib inline
from skimage.io import imread
import matplotlib.pyplot as plt
import numpy as np
bone_img = imread("ext-figures/bonegfiltslice.png")
np.random.seed(2018)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3,
figsize = (12, 4), dpi = 150)
cmap = plt.cm.nipy_spectral_r
ax1.imshow(bone_img, cmap = 'bone')
thresh_vals = [0, 70, 110, 255]
out_img = np.zeros_like(bone_img)
for i, (t_start, t_end) in enumerate(zip(thresh_vals, thresh_vals[1:])):
thresh_reg = (bone_img>t_start) & (bone_img<t_end)
ax2.hist(bone_img.ravel()[thresh_reg.ravel()],
color = cmap(i/(len(thresh_vals))))
out_img[thresh_reg] = i
ax2.set_yscale("log", nonposy='clip')
th_ax = ax3.imshow(out_img,
cmap = cmap,
vmin = 0,
vmax = len(thresh_vals))
plt.colorbar(th_ax)
<matplotlib.colorbar.Colorbar at 0x11d7415f8>
Here we end up with a goldilocks situation (Mama bear and Papa Bear), one is too low and the other is too high
fig, (ax1, ax2) = plt.subplots(1, 2, figsize = (8, 4), dpi = 300)
ax1.imshow(bone_img<thresh_vals[1], cmap = 'bone_r')
ax1.set_title('Threshold $<$ %d' % (thresh_vals[1]))
ax2.imshow(bone_img<thresh_vals[2], cmap = 'bone_r')
ax2.set_title('Threshold $<$ %d' % (thresh_vals[2]));
We can thus follow a process for ending up with a happy medium of the two
Now we apply the following steps.
%matplotlib inline
from skimage.morphology import dilation, opening, disk
from collections import OrderedDict
step_list = OrderedDict()
step_list['Strict Threshold'] = bone_img<thresh_vals[1]
step_list['Remove Small Objects'] = opening(step_list['Strict Threshold'], disk(1))
step_list['Looser Threshold'] = bone_img<thresh_vals[2]
step_list['Both Thresholds'] = 1.0*step_list['Looser Threshold'] + 1.0*step_list['Remove Small Objects']
# the tricky part keeping the between images
step_list['Connected Thresholds'] = step_list['Remove Small Objects']
for i in range(10):
step_list['Connected Thresholds'] = dilation(step_list['Connected Thresholds'] ,
disk(1.8)) & step_list['Looser Threshold']
fig, ax_steps = plt.subplots(len(step_list), 1,
figsize = (6, 4*(len(step_list))),
dpi = 150)
for i, (c_ax, (c_title, c_img)) in enumerate(zip(ax_steps.flatten(), step_list.items()),1):
c_ax.imshow(c_img, cmap = 'bone_r' if c_img.max()<=1 else 'jet')
c_ax.set_title('%d) %s' % (i, c_title))
c_ax.axis('off')
As we briefly covered last time, many measurement techniques produce quite rich data.
A pairing between spatial information (position) and some other kind of information (value). $$ \vec{x} \rightarrow \vec{f} $$
We are used to seeing images in a grid format where the position indicates the row and column in the grid and the intensity (absorption, reflection, tip deflection, etc) is shown as a different color. We take an example here of text on a page.
%matplotlib inline
from skimage.io import imread
import matplotlib.pyplot as plt
import numpy as np
from skimage.data import page
import pandas as pd
from skimage.filters import gaussian, median, threshold_triangle
page_image = page()
just_text = median(page_image, np.ones((2,2)))-255*gaussian(page_image, 20.0)
plt.imshow(page_image, cmap = 'bone');
xx, yy = np.meshgrid(np.arange(page_image.shape[1]),
np.arange(page_image.shape[0]))
page_table = pd.DataFrame(dict(x = xx.ravel(),
y = yy.ravel(),
intensity = page_image.ravel(),
is_text = just_text.ravel()>0))
page_table.sample(5)
intensity | is_text | x | y | |
---|---|---|---|---|
38169 | 174 | True | 153 | 99 |
17741 | 161 | True | 77 | 46 |
38819 | 107 | True | 35 | 101 |
24130 | 234 | True | 322 | 62 |
15369 | 120 | True | 9 | 40 |
fig, ax1 = plt.subplots(1,1)
for c_cat, c_df in page_table.groupby(['is_text']):
ax1.hist(c_df['intensity'],
np.arange(255),
label = 'Text: {}'.format(c_cat),
alpha = 0.5)
ax1.set_yscale("log", nonposy='clip')
ax1.legend();
from sklearn.metrics import roc_curve, roc_auc_score
fpr, tpr, _ = roc_curve(page_table['is_text'], page_table['intensity'])
roc_auc = roc_auc_score(page_table['is_text'], page_table['intensity'])
fig, ax = plt.subplots(1,1)
ax.plot(fpr, tpr, label='ROC curve (area = %0.2f)' % roc_auc)
ax.plot([0, 1], [0, 1], 'k--')
ax.set_xlim([0.0, 1.0])
ax.set_ylim([0.0, 1.05])
ax.set_xlabel('False Positive Rate')
ax.set_ylabel('True Positive Rate')
ax.set_title('Receiver operating characteristic example')
ax.legend(loc="lower right");
Here we can improve the results by adding information. As we discussed in the second lecture on enhancement, edge-enhancing filters can be very useful for classifying images.
def dog_filter(in_img, sig_1, sig_2):
return gaussian(page_image, sig_1) - gaussian(page_image, sig_2)
fig, (ax1, ax2) = plt.subplots(1,2, figsize = (8, 4), dpi = 150)
page_edges = dog_filter(page_image, 0.5, 10)
ax1.imshow(page_image, cmap = 'jet')
ax2.imshow(page_edges, cmap = 'jet')
page_table['edges'] = page_edges.ravel()
page_table.sample(5)
intensity | is_text | x | y | edges | |
---|---|---|---|---|---|
31063 | 234 | True | 343 | 80 | 0.101055 |
17677 | 115 | True | 13 | 46 | 0.010855 |
11181 | 110 | False | 45 | 29 | -0.041728 |
55242 | 229 | True | 330 | 143 | 0.001557 |
17891 | 218 | True | 227 | 46 | 0.053318 |
fig, ax1 = plt.subplots(1,1)
for c_cat, c_df in page_table.groupby(['is_text']):
ax1.hist(c_df['edges'],
label = 'Text: {}'.format(c_cat),
alpha = 0.5)
ax1.set_yscale("log", nonposy='clip')
ax1.legend();
from sklearn.metrics import roc_curve, roc_auc_score
fpr2, tpr2, _ = roc_curve(page_table['is_text'],
page_table['intensity']/1000.0+page_table['edges'])
roc_auc2 = roc_auc_score(page_table['is_text'],
page_table['intensity']/1000.0+page_table['edges'])
fig, ax = plt.subplots(1,1)
ax.plot(fpr, tpr, label='Intensity curve (area = %0.2f)' % roc_auc)
ax.plot(fpr2, tpr2, label='Combined curve (area = %0.2f)' % roc_auc2)
ax.plot([0, 1], [0, 1], 'k--')
ax.set_xlim([0.0, 1.0])
ax.set_ylim([0.0, 1.05])
ax.set_xlabel('False Positive Rate')
ax.set_ylabel('True Positive Rate')
ax.set_title('Receiver operating characteristic example')
ax.legend(loc="lower right");
%matplotlib inline
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
test_pts = pd.DataFrame(make_blobs(n_samples = 200, random_state = 2018)[0], columns=['x', 'y'])
plt.plot(test_pts.x, test_pts.y, 'r.')
test_pts.sample(5)
x | y | |
---|---|---|
190 | 7.504196 | -7.803367 |
96 | 8.428999 | -7.504711 |
135 | 8.877641 | -6.836999 |
73 | 5.160910 | -6.364902 |
81 | 10.318225 | -3.544903 |
from sklearn.cluster import KMeans
km = KMeans(n_clusters=2, random_state = 2018)
n_grp = km.fit_predict(test_pts)
plt.scatter(test_pts.x, test_pts.y, c = n_grp)
grp_pts = test_pts.copy()
grp_pts['group'] = n_grp
grp_pts.groupby(['group']).apply(lambda x: x.sample(5))
x | y | group | ||
---|---|---|---|---|
group | ||||
0 | 115 | 8.375861 | -6.019445 | 0 |
133 | 8.783014 | -2.409748 | 0 | |
175 | 9.376021 | -3.907714 | 0 | |
51 | 7.203155 | -7.366500 | 0 | |
40 | 9.423938 | -4.488018 | 0 | |
1 | 139 | -0.389454 | 3.579935 | 1 |
8 | -0.737558 | 2.935987 | 1 | |
123 | -1.765799 | 1.480105 | 1 | |
188 | 0.185082 | 0.862955 | 1 | |
3 | -0.350107 | 0.682015 | 1 |
We give as an initial parameter the number of groups we want to find and possible a criteria for removing groups that are too similar
What vector space to we have?
from sklearn.cluster import KMeans
km = KMeans(n_clusters = 4, random_state = 2018)
n_grp = km.fit_predict(test_pts)
plt.scatter(test_pts.x, test_pts.y, c = n_grp)
grp_pts = test_pts.copy()
grp_pts['group'] = n_grp
grp_pts.groupby(['group']).apply(lambda x: x.sample(3))
x | y | group | ||
---|---|---|---|---|
group | ||||
0 | 31 | 7.365167 | -8.297378 | 0 |
170 | 5.481257 | -7.533075 | 0 | |
71 | 7.588094 | -8.847177 | 0 | |
1 | 104 | -1.496471 | 2.038459 | 1 |
33 | -0.710215 | 2.507882 | 1 | |
119 | -1.869532 | 1.362563 | 1 | |
2 | 44 | 8.518503 | -2.729726 | 2 |
158 | 7.074534 | -4.748048 | 2 | |
7 | 8.150797 | -4.020069 | 2 | |
3 | 100 | 9.035321 | -4.372422 | 3 |
40 | 9.423938 | -4.488018 | 3 | |
163 | 9.369514 | -3.166802 | 3 |
%matplotlib inline
from skimage.io import imread
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
cortex_img = imread("ext-figures/cortex.png")[::3, ::3]/1000.0
np.random.seed(2018)
fig, (ax1) = plt.subplots(1, 1,
figsize = (8, 8), dpi = 72)
ax1.imshow(cortex_img, cmap = 'bone')
xx, yy = np.meshgrid(np.arange(cortex_img.shape[1]),
np.arange(cortex_img.shape[0]))
cortex_df = pd.DataFrame(dict(x = xx.ravel(),
y = yy.ravel(),
intensity = cortex_img.ravel()))
cortex_df.sample(5)
intensity | x | y | |
---|---|---|---|
18621 | 6.015 | 126 | 137 |
8830 | 18.961 | 55 | 65 |
23217 | 57.508 | 132 | 171 |
7605 | 5.810 | 45 | 56 |
13164 | 18.338 | 69 | 97 |
from sklearn.cluster import KMeans
km = KMeans(n_clusters = 4, random_state = 2018)
cortex_df['group'] = km.fit_predict(cortex_df[['x', 'y', 'intensity']].values)
fig, (ax1, ax2) = plt.subplots(1, 2,
figsize = (12, 8), dpi = 72)
ax1.imshow(cortex_img, cmap = 'bone')
ax2.imshow(cortex_df['group'].values.reshape(cortex_img.shape), cmap = 'gist_earth')
cortex_df.groupby(['group']).apply(lambda x: x.sample(3))
intensity | x | y | group | ||
---|---|---|---|---|---|
group | |||||
0 | 19566 | 5.843 | 126 | 144 | 0 |
18992 | 28.310 | 92 | 140 | 0 | |
21448 | 2.670 | 118 | 158 | 0 | |
1 | 10627 | 18.115 | 97 | 78 | 1 |
7125 | 26.869 | 105 | 52 | 1 | |
5103 | 0.046 | 108 | 37 | 1 | |
2 | 18766 | 17.494 | 1 | 139 | 2 |
22199 | 31.432 | 59 | 164 | 2 | |
21870 | 0.000 | 0 | 162 | 2 | |
3 | 1098 | 1.580 | 18 | 8 | 3 |
10441 | 12.067 | 46 | 77 | 3 | |
6898 | 33.534 | 13 | 51 | 3 |
Since the distance is currently calculated by $||\vec{v}_i-\vec{v}_j||$ and the values for the position is much larger than the values for the Intensity, Sobel or Gaussian they need to be rescaled so they all fit on the same axis $$\vec{v} = \left\{\frac{x}{10}, \frac{y}{10}, \textrm{Intensity}\right\}$$
km = KMeans(n_clusters = 4, random_state = 2018)
scale_cortex_df = cortex_df.copy()
scale_cortex_df.x = scale_cortex_df.x/10
scale_cortex_df.y = scale_cortex_df.y/10
scale_cortex_df['group'] = km.fit_predict(scale_cortex_df[['x', 'y', 'intensity']].values)
fig, (ax1, ax2) = plt.subplots(1, 2,
figsize = (12, 8), dpi = 72)
ax1.imshow(cortex_img, cmap = 'bone')
ax2.imshow(scale_cortex_df['group'].values.reshape(cortex_img.shape),
cmap = 'gist_earth')
scale_cortex_df.groupby(['group']).apply(lambda x: x.sample(3))
intensity | x | y | group | ||
---|---|---|---|---|---|
group | |||||
0 | 9487 | 28.010 | 3.7 | 7.0 | 0 |
22199 | 31.432 | 5.9 | 16.4 | 0 | |
3384 | 47.176 | 0.9 | 2.5 | 0 | |
1 | 14723 | 0.000 | 0.8 | 10.9 | 1 |
20250 | 0.000 | 0.0 | 15.0 | 1 | |
21877 | 0.000 | 0.7 | 16.2 | 1 | |
2 | 13381 | 26.255 | 1.6 | 9.9 | 2 |
4649 | 23.508 | 5.9 | 3.4 | 2 | |
6833 | 17.053 | 8.3 | 5.0 | 2 | |
3 | 11212 | 0.005 | 0.7 | 8.3 | 3 |
1198 | 5.308 | 11.8 | 0.8 | 3 | |
1057 | 4.233 | 11.2 | 0.7 | 3 |
km = KMeans(n_clusters = 5, random_state = 2019)
scale_cortex_df = cortex_df.copy()
scale_cortex_df.x = scale_cortex_df.x/5
scale_cortex_df.y = scale_cortex_df.y/5
scale_cortex_df['group'] = km.fit_predict(scale_cortex_df[['x', 'y', 'intensity']].values)
fig, (ax1, ax2) = plt.subplots(1, 2,
figsize = (12, 8), dpi = 72)
ax1.imshow(cortex_img, cmap = 'bone')
ax2.imshow(scale_cortex_df['group'].values.reshape(cortex_img.shape),
cmap = 'nipy_spectral')
scale_cortex_df.groupby(['group']).apply(lambda x: x.sample(3))
intensity | x | y | group | ||
---|---|---|---|---|---|
group | |||||
0 | 299 | 0.181 | 5.8 | 0.4 | 0 |
93 | 7.562 | 18.6 | 0.0 | 0 | |
1193 | 6.439 | 22.6 | 1.6 | 0 | |
1 | 20193 | 33.156 | 15.6 | 29.8 | 1 |
22053 | 31.715 | 9.6 | 32.6 | 1 | |
18599 | 34.052 | 20.8 | 27.4 | 1 | |
2 | 22564 | 14.956 | 3.8 | 33.4 | 2 |
11156 | 19.179 | 17.2 | 16.4 | 2 | |
14095 | 12.873 | 11.0 | 20.8 | 2 | |
3 | 22217 | 0.000 | 15.4 | 32.8 | 3 |
22886 | 0.000 | 14.2 | 33.8 | 3 | |
21020 | 0.000 | 19.0 | 31.0 | 3 | |
4 | 6625 | 22.726 | 2.0 | 9.8 | 4 |
3436 | 33.907 | 12.2 | 5.0 | 4 | |
2076 | 41.005 | 10.2 | 3.0 | 4 |
An approach for simplifying images by performing a clustering and forming super-pixels from groups of similar pixels.
Drastically reduced data size, serves as an initial segmentation showing spatially meaningful groups
We start with an example of shale with multiple phases
%matplotlib inline
from skimage.io import imread
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
shale_img = imread("../common/figures/shale-slice.tiff")
np.random.seed(2018)
fig, (ax1, ax2, ax3) = plt.subplots(1, 3,
figsize = (12, 4), dpi = 100)
ax1.imshow(shale_img, cmap = 'bone')
thresh_vals = np.linspace(shale_img.min(), shale_img.max(), 5+2)[:-1]
out_img = np.zeros_like(shale_img)
for i, (t_start, t_end) in enumerate(zip(thresh_vals, thresh_vals[1:])):
thresh_reg = (shale_img>t_start) & (shale_img<t_end)
ax2.hist(shale_img.ravel()[thresh_reg.ravel()])
out_img[thresh_reg] = i
ax3.imshow(out_img, cmap = 'gist_earth');
from skimage.segmentation import slic, mark_boundaries
fig, (ax1, ax2, ax3) = plt.subplots(1, 3,
figsize = (12, 4), dpi = 100)
shale_segs = slic(shale_img,
n_segments = 100,
compactness=0.1,
sigma=3.0)
ax1.imshow(shale_img, cmap = 'bone')
ax1.set_title('Original Image')
ax2.imshow(shale_segs, cmap = 'gist_earth')
ax2.set_title('Superpixels')
ax3.imshow(mark_boundaries(shale_img, shale_segs))
ax1.set_title('Superpixel Overlay')
Text(0.5,1,'Superpixel Overlay')
flat_shale_img = shale_img.copy()
for s_idx in np.unique(shale_segs.ravel()):
flat_shale_img[shale_segs==s_idx] = np.mean(flat_shale_img[shale_segs==s_idx])
fig, ax1 = plt.subplots(1,1,figsize = (8,8))
ax1.imshow(flat_shale_img, cmap = 'bone')
<matplotlib.image.AxesImage at 0x124946828>
fig, (ax1, ax2, ax3) = plt.subplots(1, 3,
figsize = (12, 4), dpi = 100)
ax1.imshow(shale_img, cmap = 'bone')
thresh_vals = np.linspace(flat_shale_img.min(), flat_shale_img.max(), 5+2)[:-1]
sp_out_img = np.zeros_like(flat_shale_img)
for i, (t_start, t_end) in enumerate(zip(thresh_vals, thresh_vals[1:])):
thresh_reg = (flat_shale_img>t_start) & (flat_shale_img<t_end)
sp_out_img[thresh_reg] = i
ax2.imshow(out_img, cmap = 'gist_earth');
ax2.set_title('Pixel Segmentation')
ax3.imshow(sp_out_img, cmap = 'gist_earth')
ax2.set_title('Superpixel Segmentation');
A more general approach is to use a probabilistic model to segmentation. We start with our image $I(\vec{x}) \forall \vec{x}\in \mathbb{R}^N$ and we classify it into two phases $\alpha$ and $\beta$
$$P(\{\vec{x} , I(\vec{x})\} | \alpha) \propto P(\alpha) + P(I(\vec{x}) | \alpha)+ P(\sum_{x^{\prime} \in \mathcal{N}} I(\vec{x^{\prime}}) | \alpha)$$Expanding on the hole filling issues examined before, a general problem in imaging is identifying regions of interest with in an image.
For samples like brains it is done to identify different regions of the brain which are responsible for different functions.
In material science it might be done to identify a portion of the sample being heated or understress.
There are a number of approaches depending on the clarity of the data and the