[2]
%load_ext watermark
%watermark -a "Romell D.Z." -u -d -p numpy,pandas,matplotlib,keras
The watermark extension is already loaded. To reload it, use: %reload_ext watermark Romell D.Z. last updated: 2019-02-22 numpy 1.16.1 pandas 0.23.4 matplotlib 2.2.2 keras 2.2.4

6.Biomedical Analysis

[1]
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import imageio
[2]
im = imageio.imread('snapshot/hand.png').astype('float64')
[3]
print('Data type:', im.dtype)
print('Min. value:', im.min())
print('Max value:', im.max())

# Plot the grayscale image
plt.imshow(im, vmin=0, vmax=255, cmap='gray')
plt.colorbar()
plt.axis('off')
plt.tight_layout()
plt.show()
Data type: float64 Min. value: 3.0 Max value: 224.0
[4]
import scipy.ndimage as ndi
hist = ndi.histogram(im, min=0, max=255, bins=256)
cdf = hist.cumsum() / hist.sum()
fig, axes = plt.subplots(2,1, sharex=True)
axes[0].plot(hist, label='Histogram')
axes[0].axvline(35,c="g")
axes[1].plot(cdf, label='CDF')
axes[1].axvline(35,c="g")
# plt.axis('off')
plt.tight_layout()
plt.show()
[5]
print(im.shape)
(480, 320)
[6]
# Create skin and bone masks
mask_bone = im >= 75
mask_skin = (im >= 45) & (im < 75)

# Plot the masks
fig, axes = plt.subplots(1,2)
axes[0].imshow(mask_skin, cmap='gray')
axes[1].imshow(mask_bone, cmap='gray')
plt.tight_layout()
plt.show()
[7]
mask_bone = im >= 70
im_bone = np.where(mask_bone, im, 0)
hist = ndi.histogram(im_bone, min=1, max=255, bins=255)
cdf = hist.cumsum() / hist.sum()
fig, axes = plt.subplots(3,1,figsize=(8,8))
axes[0].imshow(im_bone, cmap='gray')
axes[1].plot(hist)
axes[2].plot(cdf)
plt.tight_layout()
plt.show()
[8]
mask_bone = im >= 70
mask_dilate = ndi.binary_dilation(mask_bone, iterations=5)
mask_closed = ndi.binary_closing(mask_bone, iterations=5)

fig, axes = plt.subplots(1,3)
axes[0].imshow(mask_bone, cmap='gray')
axes[1].imshow(mask_dilate, cmap='gray')
axes[2].imshow(mask_closed, cmap='gray')
plt.tight_layout()
plt.show()
[9]
weights = [[0.11, 0.11, 0.11],
           [0.11, 0.11, 0.11], 
           [0.11, 0.11, 0.11]]

im_filt = ndi.convolve(im, weights)

print(im.dtype,im.size,im_filt.size)

fig, axes = plt.subplots(1,2)
axes[0].imshow(im, cmap='gray')
axes[1].imshow(im_filt, cmap='gray')
plt.tight_layout()
plt.show()
float64 153600 153600
[10]
im_s1 = ndi.gaussian_filter(im, sigma=1)
im_s3 = ndi.gaussian_filter(im, sigma=3)

fig, axes = plt.subplots(1,3)
axes[0].imshow(im >= 75, cmap='gray')
axes[1].imshow(im_s1 >= 75, cmap='gray')
axes[2].imshow(im_s3 >= 75, cmap='gray')
plt.tight_layout()
plt.show()
[11]
weights = [[1, 0, -1],
           [1, 0, -1], 
           [1, 0 ,-1]]

edges = ndi.convolve(im, weights)

plt.imshow(edges, cmap='seismic', vmin=-150, vmax=150)
plt.colorbar();
[12]
weights = [[1, 0, 1],
           [0, 0, 0], 
           [-1, -1 ,-1]]
edges = ndi.convolve(im, weights)
plt.imshow(edges, cmap='seismic', vmin=-150, vmax=150)
plt.colorbar();
[13]
weights = [[1, 0, -1],
           [1, 0, -1], 
           [1, 0 ,-1]]
bone_edges = ndi.convolve(im_bone, weights)
plt.imshow(bone_edges, cmap='seismic', vmin=-150, vmax=150)
plt.colorbar();
[14]
weights = [[1, 0, 1],
           [0, 0, 0], 
           [-1, -1 ,-1]]
edges = ndi.convolve(im_bone, weights)
plt.imshow(edges, cmap='seismic', vmin=-150, vmax=150)
plt.colorbar();
[15]
sobel_ax0 = ndi.sobel(im.astype('float64'), axis=0)
sobel_ax1 = ndi.sobel(im.astype('float64'), axis=1)

edges = np.sqrt(np.square(sobel_ax0) + np.square(sobel_ax1))

plt.imshow(edges, cmap='gray', vmax=75);
[16]
sobel_ax0 = ndi.sobel(im.astype('uint64'), axis=0)
sobel_ax1 = ndi.sobel(im.astype('uint64'), axis=1)

edges = np.sqrt(np.square(sobel_ax0) + np.square(sobel_ax1))

plt.imshow(edges, cmap='gray', vmax=75);
[17]
im_sobel_0=ndi.sobel(im_bone.astype('uint64'), axis=0)
plt.imshow(im_sobel_0, cmap='gray',vmax=75);
[18]
im_sobel_1=ndi.sobel(im_bone.astype('uint64'), axis=1)
plt.imshow(im_sobel_1, cmap='gray', vmax=75);
[19]
im_filt.min(),im_filt.max()
(6.6000000000000005, 195.69000000000003)
[20]
im_filt = ndi.median_filter(im, size=3)
range_filter = np.arange(im_filt.min(),im_filt.max(),20)
range_filter = range_filter[1:-2] # drop first & last element
print(len(range_filter))
range_filter
8
array([ 25.,  45.,  65.,  85., 105., 125., 145., 165.])
[21]
fig, axe = plt.subplots(ncols=4,nrows=2,figsize=(15,8))
labels_array = np.empty_like(range_filter)
for i, axe in enumerate(axe.flatten()):
    mask_start = np.where(im_filt > range_filter[i], 1, 0)
    mask = ndi.binary_closing(mask_start)
    labels, nlabels = ndi.label(mask)
    labels_array[i] = nlabels
    overlay = np.where(labels > 0, labels, np.nan)
    axe.imshow(im)
    axe.imshow(overlay, cmap='rainbow', alpha=0.75);
    axe.set_title('Min. Filter:%d Num. Labels: %d'% (range_filter[i],nlabels))
    
plt.tight_layout()
plt.show()
[22]
print(labels_array)
select_size = 2
select_range = labels_array.argsort()[-select_size-1:-1][::-1]
select_range
[167. 7. 24. 49. 26. 10. 4. 5.]
array([3, 4])
[23]
range_filter_2 = range_filter[select_range]
range_filter_2
array([ 85., 105.])
[24]
range_filter_2 = np.arange(range_filter_2.min(),range_filter_2.max(),2)
range_filter_2 = range_filter_2
print(len(range_filter_2))
range_filter_2
10
array([ 85.,  87.,  89.,  91.,  93.,  95.,  97.,  99., 101., 103.])
[25]
fig, axe = plt.subplots(ncols=5,nrows=2,figsize=(15,8))

im_filt = ndi.median_filter(im, size=3)    
for i, axe in enumerate(axe.flatten()):
    mask_start = np.where(im_filt > range_filter_2[i], 1, 0)
    mask = ndi.binary_closing(mask_start)
    labels, nlabels = ndi.label(mask)
    overlay = np.where(labels > 0, labels, np.nan)
    axe.imshow(im)
#     axe.axis('off')
    axe.imshow(overlay, cmap='rainbow', alpha=0.75);
    axe.set_title('Min. Filter:%d Num. Labels: %d'% (range_filter_2[i],nlabels))
    
plt.tight_layout()
plt.savefig('snapshot/biomedical_bones_detection.jpg',bbox_inches='tight')
plt.show()