Vincent GODARD - V1- 31/01/2024
Freely inspired by :
Creating Histograms https://scikit-image.org/docs/stable/auto_examples/color_exposure/plot_equalize.html
RGB to grayscale https://scikit-image.org/docs/stable/auto_examples/color_exposure/plot_rgb_to_gray.html
Sources :
Landsat Thematic Mapper (TM) data from September 10, 1987, West of Worcester, Massachusetts (from the https://clarklabs.org/download/)
scikit-image https://scikit-image.org/
matplotlib https://matplotlib.org/
Compressed file to download => here.
Upload the following file from your storage area to the data directory (TD15/data) under this notebook is located (.ipynb):
○ "h87tm4_grey.jpg";
○ "HOW87TM3.png".
This noteboock and the "data" directory will be at the same level in the tree, both under TD15.
Some libraries are already preinstalled (e.g. Matplotlib). The Skimage library is not, and must be installed.
You then need to load them (Import function).
# To install Skimage, delete the # (before the "!").
# !pip install scikit-image
# !pip install numpy
# !pip install matplotlib
## Loading libraries
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import skimage # Bibliothèque de traitement d'images
from skimage import data, img_as_float # Convert an image to floating point format.
from skimage import exposure
## Reading the TM3 image of How Hill (MA) in grayscale
HOW87TM3 = skimage.io.imread(fname='data/HOW87TM3.png', as_gray=True)
## Display image
fig, ax = plt.subplots()
plt.imshow(HOW87TM3, cmap='gray')
plt.show()
Why is it so dark?
Have a look at the histogram
## Choose your image
img = HOW87TM3
## First, let's adapt the Matplotlib settings.
plt.rcParams['font.size'] = 8 # An instance of RcParams for handling default Matplotlib values.
## Next, let's adapt the settings from scikit-image (https://scikit-image.org/docs/stable/auto_examples/color_exposure/plot_equalize.html#id4)
def plot_img_and_hist(image, axes, bins=256): # function name and parameters
"""Here's another way to leave a comment! Between two sets of "quotation marks"
Plot an image along with its histogram and cumulative histogram.
"""
# Definition of parametres
image = img_as_float(image)
ax_img, ax_hist = axes
ax_cdf = ax_hist.twinx()
# Display image
ax_img.imshow(image, cmap=plt.cm.gray)
ax_img.set_axis_off()
# Display histogram
ax_hist.hist(image.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(image, bins) # Return cumulative distribution function (cdf) for the given image.
ax_cdf.plot(bins, img_cdf, 'r')
ax_cdf.set_yticks([])
return ax_img, ax_hist, ax_cdf # special statement that you can use inside a function to send the function's result back. The return keyword is followed by some optional return values.
## Contrast stretching (Linear Stretches with Saturation)
p2, p98 = np.percentile(img, (2, 98)) # Linear Stretches with Saturation at 4%
img_rescale = exposure.rescale_intensity(img, in_range=(p2, p98))
## Equalization (assign the same number of pixels to each data level in the output image)
img_eq = exposure.equalize_hist(img)
## Adaptive Equalization [Contrast Limited Adaptive Histogram Equalization (CLAHE).]
## Algorithm for local contrast enhancement
img_adapteq = exposure.equalize_adapthist(img, clip_limit=0.03) # Clipping limit, normalized between 0 and 1 (higher values give more contrast).
## Display results
fig = plt.figure(figsize=(16, 10)) # Figure size in inches (default). Change values if necessary (columns, lines).
axes = np.zeros((2, 4), dtype=object)
axes[0, 0] = fig.add_subplot(2, 4, 1)
for i in range(1, 4):
axes[0, i] = fig.add_subplot(2, 4, 1+i, sharex=axes[0,0], sharey=axes[0,0])
for i in range(0, 4):
axes[1, i] = fig.add_subplot(2, 4, 5+i)
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(img_rescale, axes[:, 1])
ax_img.set_title('Contrast stretching')
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_eq, axes[:, 2])
ax_img.set_title('Histogram equalization')
ax_img, ax_hist, ax_cdf = plot_img_and_hist(img_adapteq, axes[:, 3])
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()
plt.show()
## Best stretch display with raw original data.
# Parameters
original = HOW87TM3
equalisation = img_adapteq # best stretch for you
fig, axes = plt.subplots(1, 2, figsize=(16, 8)) # Figure size in inches (default). Change values if necessary (columns, lines).
ax = axes.ravel()
ax[0].imshow(HOW87TM3, cmap=plt.cm.gray)
ax[0].set_title("HOW87TM3")
ax[1].imshow(img_adapteq, cmap=plt.cm.gray)
ax[1].set_title("img_adapteq")
fig.tight_layout()
plt.show()
