TP: Introduction to image processing - Radiometry¶
1. Introduction¶
For the following exercices, you need to have an install of Python 3 with numpy and matplotlib.
You need first to download the file TP.zip which contains the images you need for the session, and to put these images in the same directory as your python file.
import matplotlib.pyplot as plt
%matplotlib inline
imrgb =plt.imread("../im/parrot.png")
#the image is of type int, we cast it in float in order to do computations
imrgb=imrgb.astype(float)
# extract the three (red,green,blue) channels from y
imred = imrgb[:,:,0]
imgreen = imrgb[:,:,1]
imblue = imrgb[:,:,2]
#image size
[nrow,ncol,nch]=imrgb.shape
print(nrow,ncol,nch)
plt.imshow(imrgb)
Image display¶
#we display an image thanks to the function imshow of the pyplot library of matplotlib
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 10))
#we display the images
axes[0, 0].imshow(imrgb)
axes[0,0].set_title('Original image')
axes[0, 1].imshow(imred, cmap="gray")
axes[0,1].set_title('red channel')
axes[1, 0].imshow(imgreen, cmap="gray")
axes[1,0].set_title('green channel')
axes[1, 1].imshow(imblue, cmap="gray")
axes[1,1].set_title('blue channel')
fig.tight_layout()
It might be useful to convert the color image to gray level. This can be done by averaging the three channels, or by computing another well chosen linear combination of the coordinates R, G and B. The rgb2gray function in the image processing toolbox computes 0.2989∗R+0.5870∗G+0.1140∗B.
#we need the numpy library
import numpy as np
imgray = np.sum(imrgb,2)/3
plt.figure(figsize=(5, 5))
plt.imshow(imgray,cmap='gray')
Be careful with the function imshow() of matplotlib. For the function to work properly, it should be given an image in float/double between 0 and 1, OR an image of integers between 0 and 255.
Histograms and contrast enhancement¶
Computing histograms¶
In the following, we compute and display the gray level histogram and the cumulative histogram of an image.
The cumulative histogram of a discrete image $u$ is an increasing function defined on $\mathbb{R}$ by $$H_u(\lambda)=\frac{1}{|\Omega|}\#{\{\textbf{x};\;u(\textbf{x})\leq \lambda\}}.$$ The histogram of $u$ is the derivative of $H_u$ in the sense of distributions.
# compute histogram and cumulative histogram
imhisto, bins = np.histogram(imgray, range=(0,1), bins = 256)
imhistocum = np.cumsum(imhisto)
#we display the image, histogram and cumulative histogram
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(12, 3))
axes[0].imshow(imgray,cmap='gray')
axes[0].set_title('Parrot Image')
axes[1].bar(np.arange(0,256),imhisto)
axes[1].set_title('histogram')
axes[2].bar(np.arange(0,256),imhistocum)
axes[2].set_title('cumulative histogram')
fig.tight_layout()
Histogram equalization¶
If $u$ is a discrete image and $h_u$ its gray level distribution, histogram equalization consists in applying a contrast change $g$ (increasing function) to $u$ such that $h_{g(u)}$ is as close as possible to a constant distribution. We can compute directly $$H_u(u)*255$$.
#Image equalization
imeq = imhistocum[np.uint8(imgray*255)]
imeq=imeq.reshape(nrow,ncol)
#Display the result
plt.imshow(imeq,cmap = 'gray')
plt.show()
Apply the histogram equalization to the images ~parrot_bright.bmp~ and =parrot_dark.bmp= and plot the corresponding histograms and cumulative histograms. Comment the results and explain the observed differences.
imrgb1 = plt.imread('../im/parrot_bright.png')
imrgb2 = plt.imread('../im/parrot_dark.png')
imgray1 = imrgb1[:,:,0]
imgray2 = imrgb2[:,:,0]
[nrow,ncol] = imgray1.shape
imhisto1,bins= np.histogram(imgray1, range=(0,1), bins = 256)
imhistocum1 = np.cumsum(imhisto1)
imhisto2,bins= np.histogram(imgray2, range=(0,1), bins = 256)
imhistocum2 = np.cumsum(imhisto2)
imeq1 = imhistocum1[np.uint8(imgray1*255)]
imeq1=imeq1.reshape(nrow,ncol)
imeq2 = imhistocum2[np.uint8(imgray2*255)]
imeq2=imeq2.reshape(nrow,ncol)
#Display images
fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(12, 5))
axes[0].set_title('Bright image')
axes[0].imshow(imgrb1)
axes[1].set_title('Bright image equalized')
axes[1].imshow(imeq1,cmap = 'gray')
axes[2].set_title('Dark image')
axes[2].imshow(imgrb2)
axes[3].set_title('Dark image equalized')
axes[3].imshow(imeq2,cmap = 'gray')
fig.tight_layout()
plt.show()
Histogram specification¶
If $u$ is a discrete image and $h_u$ its gray level distribution, histogram specification consists in applying a contrast change $g$ (an increasing function) to $u$ such that $h_{g(u)}$ is as close as possible to a target discrete probability distribution $h_t$. Specification is particularly useful to compare two images of the same scene (in this case the target distribution is the histogram of the second image $v$).
- Histogram specification between two grey level images $u$ and $v$ can be computed easily by sorting the pixels of both images and by replacing each gray level in $u$ by the gray level of similar rank in $v$.
buenos1=np.double(plt.imread('../im/buenosaires3.png'))
buenos2=np.double(plt.imread('../im/buenosaires4.png'))
u = buenos1[:,:,0]
v = buenos2[:,:,1]
[nrowu,ncolu]=u.shape
[nrowv,ncolv]=v.shape
u_sort,index_u=np.sort(u,axis=None),np.argsort(u,axis=None)
[v_sort,index_v]=np.sort(v,axis=None),np.argsort(v,axis=None)
uspecifv= np.zeros(nrowu*ncolu)
uspecifv[index_u] = v_sort
uspecifv = uspecifv.reshape(nrowu,ncolu)
#Display images
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 5))
axes[0].set_title('image u')
axes[0].imshow(u,'gray')
axes[1].set_title('image v')
axes[1].imshow(v,'gray')
axes[2].set_title('image specification')
axes[2].imshow(uspecifv,'gray')
fig.tight_layout()
plt.show()
Translate the grey levels of $u$ such that it has the same mean grey level than $v$ and display the result. Is it similar to the specification of $u$ on $v$ ?
Same question by applying an affine transform to $u$ so that its mean and variance match the ones of $v$.
Midway histogram¶
The Midway histogram between two histograms $h_u$ et $h_v$ is defined from it cumulative function $H_{midway}$ : $$H_{midway}=\left( \frac{H_u^{-1}+H_v^{-1}}{2}\right)^{-1}.$$ The goal is to modify the contrast of both images $u$ and $v$ in order to give them the same intermediary grey level distribution.
u_midway=np.zeros(len(index_u))
v_midway=np.zeros(len(index_v))
u_midway[index_u] = (u_sort + v_sort)/2
v_midway[index_v] = (u_sort + v_sort)/2
u_midway = u_midway.reshape(nrowu,ncolu)
v_midway = v_midway.reshape(nrowv,ncolv)
#Display the results
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 8))
axes[0,0].set_title('image u')
axes[0,0].imshow(u,'gray')
axes[0,1].set_title('image v')
axes[0,1].imshow(v,'gray')
axes[1,0].set_title('image u_midway')
axes[1,0].imshow(u_midway,'gray')
axes[1,1].set_title('image v_midway')
axes[1,1].imshow(v_midway,'gray')
fig.tight_layout()
plt.show()
Simple transformations¶
In this exercice, you are asked to perform simple transformations on an image and find out what happens on the corresponding histogram : thresholding, affine transformation, gamma correction.
Effect of Noise on histograms¶
In the following, we want to create different noisy versions of an image $u$ and observe how the histogram $h_u$ is transformed.
Gaussian noise Write a function adding a gaussian noise $b$ to the image $u$. An image of gaussian noise of mean $0$ and of standard deviation $\sigma$ is obtained with the command ~b=sigmarandn(nrow,ncol);~ on Matlab and ~b=sigmarand(nrow,ncol,"normal");~ on Scilab. 1) Display the noisy image and its histogram for different values of $\sigma$. 2) What do you observe ? What is the relation between the histogram of $u$ and the one of $u+b$ ?
Uniform noise Same questions with $b$ a uniform additive noise.
Impulse noise Let us recall that impulse noise destroy randomly a proportion $p$ of the pixels in $u$ and replace their values by uniform random values between $0$ and $255$. Mathematically, this can be modeled as $u_b= (1-X)u+XY$, where $X$ follows a Bernouilli law of parameter $p$ and $Y$ follows a uniform law on $\{0,\dots 255\}$.
- Write a function adding impulse noise of parameter p to an image u.
/Hint/ : you can start by using the function ~rand~ to create a table I of random numbers following the uniform law on [0,1[
and then replace randomly p% of the pixels of u by a random grey levelI = rand(size(U));Ub = 255*rand(size(U)).*(I<p/100)+(I>=p/100).*U; - Display the noisy image and its histogram for different values of p. What is the relation between the histogram of u and the one of u_b ?
Image Quantization¶
Quantization¶
Image quantization consists in reducing the set of grey levels $Y = \{ y_0,\dots y_{n-1} \}$ or colors of an image $u$ into a smaller set of quantized values $\{q_0,\dots q_{p-1}\}$ ($p < n$). This operation is useful for displaying an image $u$ on a screen that supports a smaller number of colors (this is needed with a standard screen if $u$ is coded on more than 8 bits by channel).
A quantization operator $Q$ is defined by the values $(q_i)_{i=0, \dots p-1}$ and $(t_j)_{j=0,\dots p}$ such that $$ t_0 \leq q_0 \leq t_1 \leq q_1 \leq \dots q_{p-1} \leq t_p,\text{ and } Q(\lambda)=q_i \text{ if } t_i \leq \lambda < t_{i+1}.$$
Uniform Quantization Uniform quantization consists in dividing the set $Y$ in $p$ regular intervals.
- Use uniform quantization on the grey level version of ~lena~ (try different numbers $K$ of grey levels) and display the result. For which value of $K$ do you start to see a difference with the original image ?
u = plt.imread('lena.png')
for K in range(50,1,-10):
uquantif = (np.floor(u*(K-1)+0.5))/K# quantization on K=10 grey levels
plt.imshow(uquantif,'gray')
plt.show()
Histogram-based Quantization
This consists in choosing $t_i=\min \{\lambda; H_u(\lambda) \geq \frac{i}{p} \}$, and the $q_i$ are defined as the barycenters of the intervals $[t_i,t_{i+1}].$
1) Show that this boils down to an histogram equalization followed by a uniform quantization
2) Apply this quantization on the grey level version of ~lena.bmp~ and display the result. Same question on the limit value $K$ for which we perceive a difference with the original image.
Lloyd-Max quantization This quantization consists in minimizing the least square error $$LSE((q_i)_{i=1\dots p-1},(t_i)_{i=1\dots p})= \sum_{i=0}^{p-1} \int_{t_i}^{t_{i+1}} h(\lambda) (\lambda -q_i)^2.$$ 1) Write the optimality conditions that should be satisfied by the solution $\{(\widehat{q_i}),(\widehat{t_i})\}$. 2) Write a function which minimizes the least square error by alternatively minimizing in $(q_i)_{i=0, \dots p-1}$ and $(t_j)_{j=0,\dots p}$. 3) Apply this quantization on the grey level version of ~lena.bmp~ for different values of $K$ and display the result. Comment.
Dithering¶
Dithering consists in adding intentionnally noise to an image before quantization. For instance, it can be used to convert a grey level image to black and white in such a way that the density of white dots in the new image is an increasing function of the grey level in the original image. This is particularly useful for impression or displaying.
Let us explain how dithering works in the case of 2 grey levels (binarization). All grey levels smaller than a value $\lambda$ are replaced by $0$ and those greater than $\lambda$ are replaced by $255$. If we add a i.i.d. noise $B$ of density $p_B$ to $u$ before the binarization, then at the pixel $x$ we get $$P[u(x) + B(x) > \lambda] = P[B(x) > \lambda - u(x) ] = \int_{\lambda - u(x)}^{+\infty} p_B(s)ds,$$ which is an increasing function of the value $u(x)$. The probability that $x$ turns white in the dithered image is thus an increasing function of its original grey level.
- Perform dithering to reduce the number of grey levels of ~lena.bmp~. Compare the result with a quantization without dithering.