For the following exercices, you need Python 3 with some basic librairies (see below). All images necessary for the session are available here.
If you use your own Python 3 install, you should download the images, put them in a convenient directory and update the path in the next cell.
path = '../im/'
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
A color image is made of three channels : red, green and blue. A color image in $\mathbb{R}^{N\times M}$ is stored as a $N\times M\times 3$ matrix.
Be careful with the functions plt.imread()
and plt.imshow()
of matplotlib
.
plt.imread()
reads png images as numpy arrays of floating points between 0 and 1, but it reads jpg or bmp images as numpy arrays of 8 bit integers.
In this practical session, we assume images encoded as floating point values between 0 and 1, so if you load a jpg or bmp file you must convert the image to float type and normalize its values.
If 'im' is an image encoded as a double numpy array, plt.imshow(im)
will display all values above 1 in white and all values below 0 in black. If the image 'im' is encoded on 8 bits though, plt.imshow(im)
will display 0 in black and 255 in white.
imrgb = plt.imread(path+'parrot.png')
[nrow,ncol,nch] = imrgb.shape
print(nrow,ncol,nch)
495 495 3
You can use plt.imshow()
to display the 3D numpy array imrgb
as an image:
plt.figure(figsize=(7, 7))
plt.imshow(imrgb)
plt.show()
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. First we try with simple averaging $$I_{gs}=(R+G+B)/3$$
imgray = np.mean(imrgb,2)
plt.figure(figsize=(7, 7))
plt.imshow(imgray,cmap='gray',vmin=0,vmax=1)
plt.show()
1. Now use a custom weighted averaging of the three channels, that reflects better human perception: $$I_{gs}=0.21 R + 0.72 G + 0.07 B$$
# to do: question 1
imgray2 = 0.21*imrgb[:,:,0] + 0.72*imrgb[:,:,1] + 0.07*imrgb[:,:,2]
# autre solution : imgray2 = np.average(imrgb,axis=2,weights=[0.21,0.71,0.07])
plt.figure(figsize=(7, 7))
plt.imshow(imgray2,cmap='gray',vmin=0,vmax=1)
plt.show()
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.
a. We compute the histogram of the image imgray
:
imhisto, bins = np.histogram(imgray, range=(0,1), bins = 256)
imhisto = imhisto/np.sum(imhisto)
b. We now compute the corresponding cumulative histogram thanks to the function np.cumsum()
which cumulates the values of a vector from left to right:
imhistocum = np.cumsum(imhisto)
c. We display the image, histogram and cumulative histogram:
values = (bins[1:]+bins[:-1])/2
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(15, 5))
axes[0].imshow(imgray,cmap='gray',vmin=0,vmax=1)
axes[0].set_title('parrot image, gray level')
axes[1].bar(values,imhisto,width=1/256)
axes[1].set_title('histogram')
axes[2].bar(values,imhistocum,width=1/256)
axes[2].set_title('cumulative histogram')
fig.tight_layout()
plt.show()
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.$$
To this aim, we can apply directly the vector imhistocum
(which can be seen as a function from $\{0,\dots,255\}$ into $[0,1]$) to the numpy array imgray
. Since imgray
has values between $0$ and $1$, it is necessary to multiply it by $255$ and cast it as a 8-bit array:
imgray_8bits = np.uint8(255*imgray) # changement de l'encodage des valeurs de niveaux de gris de imgray
imeq = imhistocum[imgray_8bits]
The resulting image imeq
has the same shape as imgray
. To better understand what operation does imhistocum[imgray_8bits]
, here is an example with small size matrices :
a = np.array([[3,4,5],[3,1,0]]) # matrix of integers, analog to imgray_8bits
print("a=",a)
b = np.array([7,8,9,10,11,12]) # vector, analog to imhistocum
print("b=",b)
c = b[a]
print("c=",c) # c has same shape as a
# each value in a is replaced by the corresponding value : 3 becomes 10, because b[3]=10,
# 4 becomes 11 because b[4]=11, and so on.
a= [[3 4 5] [3 1 0]] b= [ 7 8 9 10 11 12] c= [[10 11 12] [10 8 7]]
We can now display the resulting equalized image:
plt.figure(figsize=(7, 7))
plt.imshow(imeq,cmap = 'gray',vmin=0,vmax=1)
plt.show()
2. Compute and plot also the corresponding histograms and cumulative histograms of the equalized image.
def comp_histos(im):
imhisto, bins = np.histogram(im, range=(0,1), bins = 256)
imhisto = imhisto/np.sum(imhisto)
imhistocum = np.cumsum(imhisto)
return imhisto, bins, imhistocum
def plot_histos(im):
imhisto, bins, imhistocum = comp_histos(im)
values = (bins[1:]+bins[:-1])/2
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(15, 5))
axes[0].imshow(im,cmap='gray',vmin=0,vmax=1)
axes[1].bar(values,imhisto,width=1/256)
axes[1].set_title('histogram')
axes[2].bar(values,imhistocum,width=1/256)
axes[2].set_title('cumulative histogram')
fig.tight_layout()
plt.show()
plot_histos(imeq)
3. Now, apply the previous histogram equalization to the two images parrot_bright.png
and parrot_dark.png
, plot the equalized images, the corresponding histograms and cumulative histograms. Comment the results and explain the observed differences.
# we first read and display 'parrot_bright.png'
imrgb_bright = plt.imread(path+'parrot_bright.png')
imgray_bright = imrgb_bright[:,:,0]
plt.imshow(imgray_bright,cmap='gray',vmin=0,vmax=1)
plt.show()
# Now we define a function to compute the equalization of an image:
def equalize_image(imgray):
imhisto, bins, imhistocum = comp_histos(imgray)
imgray_8bits = np.uint8(255*imgray) # changement de l'encodage des valeurs de niveaux de gris de imgray
imeq = imhistocum[imgray_8bits]
return imeq
# We apply this function to imgray_bright and display the results with 'plot_histos'
imeq_bright = equalize_image(imgray_bright)
plot_histos(imeq_bright)
# Now we do the same with 'parrot_dark.png'
imrgb_dark = plt.imread(path+'parrot_dark.png')
imgray_dark = imrgb_dark[:,:,0]
plt.imshow(imrgb_dark, cmap='gray',vmin=0,vmax=1)
plt.show()
imeq_dark = equalize_image(imgray_dark)
plot_histos(imeq_dark)
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$).
We start by reading our two images $u$ and $v$:
buenos1=plt.imread(path+'buenosaires4.png')
buenos2=plt.imread(path+'buenosaires3.png')
u = buenos1[:,:,0]
v = buenos2[:,:,0]
[nrowu,ncolu]=u.shape
[nrowv,ncolv]=v.shape
Now, 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$:
index_u = np.argsort(u,axis=None)
u_sort = np.sort(u,axis=None)
index_v = np.argsort(v,axis=None)
v_sort = np.sort(v,axis=None)
uspecifv = np.zeros(nrowu*ncolu)
uspecifv[index_u] = v_sort
uspecifv = uspecifv.reshape(nrowu,ncolu)
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(15, 5))
axes[0].set_title('image u')
axes[0].imshow(u,'gray',vmin=0,vmax=1)
axes[1].set_title('image v')
axes[1].imshow(v,'gray',vmin=0,vmax=1)
axes[2].set_title('image specification')
axes[2].imshow(uspecifv,'gray',vmin=0,vmax=1)
fig.tight_layout()
plt.show()
4. Try to 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$ ?
# First we compute the "translation value", i.e. the difference between the mean grey levels of the two images:
diff_mean_uv = np.mean(v-u)
# Now we apply the translation:
u_trans = u + diff_mean_uv
# The operation may have produced values outside the range [0,1], so we need to do a correction:
u_trans = np.maximum(0,np.minimum(1,u_trans))
# We plot the resulting image, together with the specification of u on v:
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(15, 5))
axes[0].set_title('translated image')
axes[0].imshow(np.tile(u_trans[:,:,None],(1,1,3)),'gray',vmin=0,vmax=1)
axes[1].set_title('specification')
axes[1].imshow(uspecifv,'gray',vmin=0,vmax=1)
fig.tight_layout()
plt.show()
5. Same question by applying an affine transform to $u$ so that its mean and variance match the ones of $v$.
# We can first multiply u by the quotient of standard deviations between v and u, in order to get equal variances:
u_aff = u * (np.std(v)/np.std(u))
# Now we can do the same translation as previously:
u_aff = u_aff + np.mean(v-u_aff)
# Correction to fit in the range [0,1]:
u_aff = np.maximum(0,np.minimum(1,u_aff))
# We plot the resulting image, together with the specification of u on v:
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(15, 5))
axes[0].set_title('affine transformed image')
axes[0].imshow(u_aff,'gray',vmin=0,vmax=1)
axes[1].set_title('specification')
axes[1].imshow(uspecifv,'gray',vmin=0,vmax=1)
fig.tight_layout()
plt.show()
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=(14, 10))
axes[0,0].set_title('image u')
axes[0,0].imshow(u,'gray',vmin=0,vmax=1)
axes[0,1].set_title('image v')
axes[0,1].imshow(v,'gray',vmin=0,vmax=1)
axes[1,0].set_title('image u_midway')
axes[1,0].imshow(u_midway,'gray',vmin=0,vmax=1)
axes[1,1].set_title('image v_midway')
axes[1,1].imshow(v_midway,'gray',vmin=0,vmax=1)
fig.tight_layout()
plt.show()
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.
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
6 a) 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
sigma*np.random.randn(nrow,ncol)
6 b) Display the noisy image and its histogram for different values of $\sigma$.
6 c) What do you observe ? What is the relation between the histogram of $u$ and the one of $u+b$ ?
imgray = np.mean(plt.imread(path+'parrot.png'),2)
def add_gaussian_noise(u,sigma):
nrow, ncol = u.shape
res = u + sigma * np.random.randn(nrow,ncol)
res = np.maximum(0,np.minimum(1,res))
return res
for sigma in [0,0.05,0.1,0.2,0.5]:
imgray_noise = add_gaussian_noise(imgray,sigma)
plot_histos(imgray_noise)
Uniform noise 7. Same questions with $b$ a uniform additive noise.
imgray = np.mean(plt.imread(path+'parrot.png'),2)
def add_uniform_noise(u,sigma):
alpha = np.sqrt(3)*sigma
nrow, ncol = u.shape
res = u + (np.random.rand(nrow,ncol)-.5)*2*alpha
res = np.maximum(0,np.minimum(1,res))
return res
for sigma in [0,0.05,0.1,0.2,0.5]:
imgray_noise = add_uniform_noise(imgray,sigma)
plot_histos(imgray_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\}$.
8 a) 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 tab
of random numbers following the uniform law on $[0,1[$
tab = np.random.rand(u.shape[0],u.shape[1])
and then replace randomly $p\%$ of the pixels of $u$ by a random grey level
ub = 255*np.random.rand(u.shape[0],u.shape[1])*(tab<p/100)+(tab>=p/100)*u;
8 b) 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 ?
def impulse_noise(u,p):
nrow, ncol = u.shape
Y = np.random.rand(nrow,ncol)
tab = np.random.rand(nrow,ncol)
X = np.zeros((nrow,ncol))
X[tab<p] = 1
res = (1-X)*u+X*Y
return res
imgray = np.mean(plt.imread(path+'parrot.png'),2)
for p in [0,.1,.25,.5,.75]:
imgray_noise = impulse_noise(imgray,p)
plot_histos(imgray_noise)
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.
u = np.mean(plt.imread(path+'simpson512.png'),2)
K = 20
# 1e méthode
u_quant = np.zeros(u.shape)
for k in range(K):
test = (u>k/K) & (u<=(k+1)/K)
u_quant[test] = (k+.5)/K
# 2e méthode
u_quant = (np.floor(K*u)+.5)/K
# plusieurs valeurs de K
tabK = [50,30,20,10,5,2]
nK = len(tabK)
plt.figure(figsize=(16, 10))
for i, K in enumerate(tabK):
u_quant = (np.floor(K*u)+.5)/K
plt.subplot(2,np.ceil(nK/2),i+1)
plt.imshow(u_quant,cmap='gray',vmin=0,vmax=1)
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}].$
10 a) Show that this boils down to an histogram equalization followed by a uniform quantization
10 b) Apply this quantization on a gray level image and display the result. Same question on the limit value $K$ for which we perceive a difference with the original image.
# correction
# we define a function to perform the quantization of an image u given the values of t_i and q_i :
def quantize(u,t,q):
K = len(q)
u_quant = np.zeros(u.shape)
for k in range(K):
test = (u>=t[k]) & (u<t[k+1])
u_quant[test] = q[k]
return u_quant
def histogram_quantize(u,K):
nrow, ncol = u.shape
# we need to compute t_i = H_u^-1(i/K)
# First we sort the values of u:
S_u = np.sort(u,axis=None)
# S_u gives the sorted values of u. We have S_u[j]=lambda <=> j/|Omega| = H_u(lambda)
# So we have : t_i = H_u^-1(i/K) = S_u[|Omega|*i/K]
# So first we need to compute the indices j = |Omega|*i/K for i=0,1,...,K :
tab_i = np.arange(K+1)/K
tab_j = ((nrow*ncol-1)*tab_i).astype(int)
# now we compute the t_i = S_u[j]
t = S_u[tab_j]
# finally we compute the q_i = (t_{i+1}+t_i)/2 :
q = .5 * (t[1:]+t[:-1])
# finally we compute the quantization given t and q :
u_quant = quantize(u,t,q)
return u_quant
imgray = np.mean(plt.imread(path+'simpson512.png'),2)
plt.figure(figsize=(16, 10))
tabK = [50,30,20,10,5,2]
nK = len(tabK)
for i, K in enumerate(tabK):
im_quant = histogram_quantize(imgray,K)
plt.subplot(2,np.ceil(nK/2),i+1)
plt.title('K='+str(K))
plt.imshow(im_quant,cmap='gray',vmin=0,vmax=1)
plt.tight_layout()
plt.show()
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.$$
It is equivalent to the algorithm Kmeans
in one dimension.
11 a) Write the optimality conditions that should be satisfied by the solution $\{(\widehat{q_i}),(\widehat{t_i})\}$.
11 b) 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}$.
11 c) Apply this quantization on the previous gray level image for different values of $K$ and display the result. Comment.
N.B. We assume that $t_0=0$ and $t_p=1$ are fixed.
The optimality conditions are :
$$ \frac{\partial LSE}{\partial q_j}=0
\qquad\Leftrightarrow\qquad q_j=\frac{\int_{t_j}^{t_{j+1}}\lambda h(\lambda)d\lambda}{\int_{t_j}^{t_{j+1}}h(\lambda)d\lambda}\quad\forall j=0,1,\ldots,p-1$$
and
$$ \frac{\partial LSE}{\partial t_j}=0
\qquad\Leftrightarrow\qquad t_j=\frac{q_j+q_{j+1}}2\quad\forall j=1,\ldots,p-1$$
We can iterate on these two formulas in order to optimize alternatively over $q_j$ and $t_j$. Note that we can use the cumulative histogram to compute one of the two integrals, since $\int_{t_j}^{t_{j+1}}h(\lambda)d\lambda = H_u(t_{j+1})-H_u(t_j)$. Similarly for the other integral, we can pre-compute the corresponding cumulative sum (denoted imcumesp
in the code) so that we can just compute it as a difference during the iterations.
def Lloyd_Max_quantize(u,K,plot=False):
imhisto, bins, imhistocum = comp_histos(u)
imcumesp = np.cumsum(np.arange(256)*imhisto/256)
t = np.arange(K+1)/K
test = True
eps = 1e-6
k = 0
while test:
k+=1
indt = (255*t).astype(np.int)
q = np.diff(imcumesp[indt]) / np.diff(imhistocum[indt])
t_old = t.copy()
t[1:-1] = .5 * (q[1:]+q[:-1])
test = (np.linalg.norm(t-t_old) > eps)
if K==2:
# N.B. we prefer to use values 0 and 1 for a binary image
q[0], q[1] = 0, 1
u_quant = quantize(u,t,q)
if plot:
plt.figure(figsize=(16, 7))
plt.subplot(1,2,1)
plt.title('Lloyd-Max quantization with ' + str(K) + ' levels')
plt.imshow(u_quant,cmap='gray',vmin=0,vmax=1)
plt.subplot(1,2,2)
plt.title('histogram of u, points t_i (in blue) and q_i (in red)')
plt.plot(np.arange(256)/256,imhisto)
plt.plot(t,np.zeros(t.shape),'.b')
plt.plot(q,np.zeros(q.shape),'+r')
plt.show()
return u_quant
imgray = np.mean(plt.imread(path+'simpson512.png'),2)
u_quant = Lloyd_Max_quantize(imgray,K=10,plot=True)
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.
12 a) Perform dithering in order to quantize a gray level image on 10 levels (you can add a small Gaussian noise of std 5/255 for instance). Compare the result with the previous quantizations without dithering.
12 b) Try with different levels of noise.
def dither_quantize(u,K,sigma):
u = u + sigma * np.random.randn(u.shape[0],u.shape[1])
u_quant = Lloyd_Max_quantize(u,K)
return u_quant
def demo_dither(u,K,sigma):
im_quant = Lloyd_Max_quantize(u,K)
plt.figure(figsize=(14, 7))
plt.subplot(1,2,1)
plt.title('Lloyd-Max quantization with '+str(K)+' levels')
plt.imshow(im_quant,cmap='gray',vmin=0,vmax=1)
im_dither_quant = dither_quantize(u,K,sigma)
plt.subplot(1,2,2)
plt.title('same, with dithering (Gaussian noise of std {:f})'.format(sigma))
plt.imshow(im_dither_quant,cmap='gray',vmin=0,vmax=1)
plt.show()
imgray = np.mean(plt.imread(path+'simpson512.png'),2)
demo_dither(imgray,K=10,sigma=5/255)
imgray = np.mean(plt.imread(path+'lighthouse.bmp'),2) / 255
demo_dither(imgray,K=2,sigma=10/255)