In [1]:
# Author: Qinghui Liu, for my master thesis project,
# Date: 2017-03
# basic image processing and analysis by python
# ####1. Spatial Filters - linear filters and non-linear filters ####
# linear ones include mean, laplacian, and laplacian of gaussian
# non-linear ones include median, maximum, minimum, sobel, prewitt and canny
# Four padding approaches: Zero padding, constant, nearest neighbor and reflect paddings.
# ---mean filter
import numpy as np
import skimage
from skimage import feature
import skimage.io as sio
import scipy.misc
import scipy.ndimage as sn
from scipy.misc.pilutil import Image
from matplotlib import pyplot
#import matplotlib.pyplot as plt
# plotting inline in Jupter Notebook
%matplotlib inline
#matplotlib.rcParams['font.size'] = 8
# open image and convert it to grayscale
img_dir = '/Users/liuqh/Desktop/keras/data/train/polyp/_0_1452.jpeg'
#a = Image.open(img_dir).convert('L')
a = sio.imread(img_dir)
a = skimage.color.rgb2gray(a)
pyplot.subplot(1,2,1)
pyplot.title('input img')
pyplot.imshow(a,cmap='gray')
# initializeing the filter of size 5x5
# divided by 25 for normalization
k = np.ones((5,5))/25
# perform convolution
b = sn.filters.convolve(a,k)
pyplot.subplot(1,2,2)
pyplot.imshow(b,cmap='gray') # try cmap = 'bone_r' or other parameters
pyplot.title('mean filter output')
pyplot.show()
# convert ndarray to an image
b = scipy.misc.toimage(b)
b.save('mean_output.png')
In [2]:
### Median Filter
# median filter - one popular non-linear filter
b_median = scipy.ndimage.filters.median_filter(a, size=5, footprint=None, output=None,
mode ='reflect', cval=0.0, origin=0)
b_median = scipy.misc.toimage(b_median)
pyplot.subplot(1,2,1)
pyplot.imshow(b_median,'gray')
pyplot.title('median filter output')
b_median.save('b_median.png')
### Max Filter
# this filter enhances the bright points
b_max = sn.filters.maximum_filter(a, size=5,
footprint=None,
output=None,
mode ='reflect',
cval=0.0, origin=0)
b_max = scipy.misc.toimage(b_max)
pyplot.subplot(1,2,2)
pyplot.imshow(b_max,'gray')
pyplot.title('max filter output')
b_max.save('b_max.png')
In [3]:
### Min Filter
# this filter enhances the darkest points
b_min = sn.filters.minimum_filter(a, size=5,
footprint=None,
output=None,
mode ='reflect',
cval=0.0, origin=0)
b_min = scipy.misc.toimage(b_min)
pyplot.subplot(1,2,1)
pyplot.imshow(b_min,'gray')
pyplot.title('min filter output')
b_min.save('b_min.png')
### Edge detection
# Sobel, and Prewitt filters are used to enchance all edges
# horizontal or vertical - sobel or prewitt just enhance all vertical or horizontal endges
from skimage import filter
b_edge = filter.sobel(a) # try to use sobel_v(a) or sobel_h(a)
pyplot.subplot(1,2,2)
pyplot.imshow(b_edge,'gray')
pyplot.title('sobel filter output')
b_edge = scipy.misc.toimage(b_edge)
b_edge.save('b_edge.png')
In [4]:
# prewitt and hprewitt filters
b_prewitt = filter.prewitt(a,mask = None)
pyplot.subplot(1,2,1)
pyplot.imshow(b_prewitt,'gray')
pyplot.title('prewitt filter output')
b_prewitt = scipy.misc.toimage(b_prewitt)
b_prewitt.save('b_prewitt.png')
b_hprewitt = filter.prewitt_h(b_min,mask = None)
pyplot.subplot(1,2,2)
pyplot.imshow(b_hprewitt,'gray')
pyplot.title('hprewitt filter output')
b_hprewitt = scipy.misc.toimage(b_hprewitt)
b_hprewitt.save('b_hprewitt.png')
In [5]:
# canny and laplace filters
b_canny = feature.canny(a, sigma=0.1)
pyplot.subplot(1,2,1)
pyplot.imshow(b_canny,'gray')
pyplot.title('canny filter output')
b_canny = scipy.misc.toimage(b_canny)
b_canny.save('b_canny.png')
#b_laplace = skimage.filters.laplace(a,ksize = 3)
b_laplace = sn.filters.laplace(a,mode='reflect')
pyplot.subplot(1,2,2)
pyplot.imshow(b_laplace,'gray')
pyplot.title('laplace filter output')
b_laplace = scipy.misc.toimage(b_laplace)
b_laplace.save('b_laplace.png')
In [6]:
# Histogram Equalization
# refer to link: http://scikit-image.org/docs/dev/auto_examples
from skimage import data, img_as_float
from skimage import exposure
def plot_img_and_hist(img, axes, bins=256):
"""Plot an image along with its histogram and cumulative histogram.
"""
img = img_as_float(img)
ax_img, ax_hist = axes
ax_cdf = ax_hist.twinx()
# Display image
ax_img.imshow(img, cmap='gray')
ax_img.set_axis_off()
ax_img.set_adjustable('box-forced')
# Display histogram
ax_hist.hist(img.ravel(), bins=bins, histtype='step', color='black')
ax_hist.ticklabel_format(axis='y', style='scientific', scilimits=(0, 0))
ax_hist.set_xlabel('Pixel intensity')
ax_hist.set_xlim(0, 1)
ax_hist.set_yticks([])
# Display cumulative distribution
img_cdf, bins = exposure.cumulative_distribution(img, bins)
ax_cdf.plot(bins, img_cdf, 'r')
ax_cdf.set_yticks([])
return ax_img, ax_hist, ax_cdf
In [7]:
img = a
# Contrast stretching
p2, p98 = np.percentile(img, (2, 98))
img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))
# Equalization
img_eq = exposure.equalize_hist(img)
# Adaptive Equalization
img_adapteq = exposure.equalize_adapthist(img, clip_limit=0.03)
In [8]:
# Display results
fig = pyplot.figure(figsize=(8, 5))
axes = np.zeros((2, 3), dtype=np.object)
axes[0, 0] = fig.add_subplot(2, 3, 1)
for i in range(1, 3):
axes[0, i] = fig.add_subplot(2, 3, 1+i, sharex=axes[0,0], sharey=axes[0,0])
for i in range(0, 3):
axes[1, i] = fig.add_subplot(2, 3, 4+i)
#ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0])
#ax_img.set_title('Low contrast image')
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_rescale, axes[:, 0])
y_min, y_max = ax_hist.get_ylim()
ax_img.set_title('Contrast stretching')
ax_hist.set_ylabel('Number of pixels')
ax_hist.set_yticks(np.linspace(0, y_max, 5))
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_eq, axes[:, 1])
ax_img.set_title('Histogram equalization')
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_adapteq, axes[:, 2])
ax_img.set_title('Adaptive equalization')
ax_cdf.set_ylabel('Fraction of total intensity')
ax_cdf.set_yticks(np.linspace(0, 1, 5))
# prevent overlap of y-axis labels
fig.tight_layout()
pyplot.show()
In [9]:
# Power law transformation
# t(i,j) = kI(i,j)^r
import math, numpy
im = Image.open(img_dir).convert('L')
imp = scipy.misc.fromimage(im)
gamma = 0.2 # try to use other numbers to text 0.5, 1, 2, 5, etc
imp1 = imp.astype(float)
imp3 = numpy.max(imp1)
imp2 = imp1/imp3
# compute gamma-correction
imp3 = numpy.log(imp2)*gamma
# perform gamma-correction
c = numpy.exp(imp3)*255.0
# convert c to type int
c1 = c.astype(int)
# convert c1 from ndarray to image
im_pl = scipy.misc.toimage(c1)
im_pl.save('b_powerlaw.png')
In [10]:
## image inverse transfomation
# t(i,j) = L - 1 -I(i,j), transfrom dark intensities to bright intensities
# vice versa
im2 = 255-imp
im3 = scipy.misc.toimage(im2)
im3.save('b_invers.png')
pyplot.subplot(1,2,1)
pyplot.imshow(im_pl,'gray')
pyplot.title('img power low output')
pyplot.subplot(1,2,2)
pyplot.imshow(im3,'gray')
pyplot.title('img inverse output')
Out[10]:
In [11]:
# log transformation
# t(i,j) = k* log(1+I(i,j)) where
# k = (L-1)/log(1+|I_max|), I_max is maximum magnitude
p1=imp1.astype(float)
p2=numpy.max(p1)
c=(255.0*numpy.log(1+p1))/numpy.log(1+p2)
c1=c.astype(int)
im_log = scipy.misc.toimage(c1)
im_log.save('b_logTrans.png')
pyplot.subplot(1,2,2)
pyplot.imshow(im_log,'gray')
pyplot.title('img log trans output')
Out[11]:
In [12]:
# Gamma and log contrast adjustment
# Gamma
gamma_corrected = exposure.adjust_gamma(img, 2)
# Logarithmic
logarithmic_corrected = exposure.adjust_log(img, 1)
In [13]:
# Display results
fig = pyplot.figure(figsize=(8, 5))
axes = np.zeros((2, 3), dtype=np.object)
axes[0, 0] = pyplot.subplot(2, 3, 1, adjustable='box-forced')
axes[0, 1] = pyplot.subplot(2, 3, 2, sharex=axes[0, 0], sharey=axes[0, 0],
adjustable='box-forced')
axes[0, 2] = pyplot.subplot(2, 3, 3, sharex=axes[0, 0], sharey=axes[0, 0],
adjustable='box-forced')
axes[1, 0] = pyplot.subplot(2, 3, 4)
axes[1, 1] = pyplot.subplot(2, 3, 5)
axes[1, 2] = pyplot.subplot(2, 3, 6)
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img, axes[:, 0])
ax_img.set_title('Low contrast image')
y_min, y_max = ax_hist.get_ylim()
ax_hist.set_ylabel('Number of pixels')
ax_hist.set_yticks(np.linspace(0, y_max, 5))
ax_img, ax_hist, ax_cdf = plot_img_and_hist(gamma_corrected, axes[:, 1])
ax_img.set_title('Gamma correction')
ax_img, ax_hist, ax_cdf = plot_img_and_hist(logarithmic_corrected, axes[:, 2])
ax_img.set_title('Logarithmic correction')
ax_cdf.set_ylabel('Fraction of total intensity')
ax_cdf.set_yticks(np.linspace(0, 1, 5))
# prevent overlap of y-axis labels
fig.tight_layout()
pyplot.show()
In [24]:
# adapting gray-scale filters to RGB images
# each_channel, pass each of RGB channels to the filter
# hsv_vaule, convert RGB to HSV and pass the value channel to the filter
# the result is inserted back to HSV and then converted back to RGB
from skimage.color.adapt_rgb import adapt_rgb, each_channel, hsv_value
from skimage import filters
from skimage.exposure import rescale_intensity
@adapt_rgb(each_channel)
def sobel_each(image):
return filters.sobel(image)
@adapt_rgb(hsv_value)
def sobel_hsv(image):
return filters.sobel(image)
image = sio.imread(img_dir) # rgb imaaage
fig = pyplot.figure(figsize=(8, 4))
ax_each = fig.add_subplot(131, adjustable='box-forced')
ax_hsv = fig.add_subplot(132, sharex=ax_each, sharey=ax_each,
adjustable='box-forced')
ax_orig = fig.add_subplot(133, adjustable='box-forced')
# We use 1 - sobel_each(image)
# but this will not work if image is not normalized
ax_each.imshow(rescale_intensity(1 - sobel_each(image)))
ax_each.set_xticks([]), ax_each.set_yticks([])
ax_each.set_title("Sobel filter \n on each RGB chan")
# We use 1 - sobel_hsv(image) but this will not work if image is not normalized
ax_hsv.imshow(rescale_intensity(1 - sobel_hsv(image)))
ax_hsv.set_xticks([]), ax_hsv.set_yticks([])
ax_hsv.set_title("Sobel filter \n on each HSV ")
ax_orig.imshow(image)
ax_orig.set_xticks([]), ax_orig.set_yticks([])
ax_orig.set_title("Original RGB image")
Out[24]:
In [25]:
## Thresholding
# create a binary image from a grayscale image
# http://en.wikipedia.org/wiki/Otsu’s_method
from skimage.filters import threshold_otsu,threshold_isodata, threshold_li
#image = skimage.color.rgb2gray(image)
#image= img_eq#.astype(float)
thresh = threshold_otsu(image)
#thresh = threshold_li(image)
#thresh = threshold_isodata(image)
binary = image < thresh # try to test other thresholds
fig, axes = pyplot.subplots(ncols=3,figsize=(8, 2.5))
ax = axes.ravel()
ax[0] = pyplot.subplot(1, 3, 1, adjustable='box-forced')
ax[1] = pyplot.subplot(1, 3, 2)
ax[2] = pyplot.subplot(1, 3, 3, sharex=ax[0], sharey=ax[0], adjustable='box-forced')
ax[0].imshow(image, cmap=pyplot.cm.gray)
ax[0].set_title('Original')
ax[0].axis('off')
ax[1].hist(image.ravel(), bins=256)
ax[1].set_title('Histogram')
ax[1].axvline(thresh, color='r')
# binary.astype(np.uint8), np.float, np.uint16
ax[2].imshow(binary.astype(np.float), cmap=pyplot.cm.gray)
ax[2].set_title('Thresholded')
ax[2].axis('off')
pyplot.show()
In [ ]: