1. notebooks
  2. python

Introduction to image processing - Radiometry

Introduction

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.

For some parts of the session (cells with commands written as todo_something...), you are supposed to code by yourself.

In [1]:
path = '../im/'
In [2]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

The following line will be used to import the solutions of the practical session. Do not use it for the moment.

In [3]:
from TP_radiometrie import *

Load and display a color image

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 are 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 to $[0,1]$.

  • 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.</span>

In [4]:
imrgb =plt.imread(path+"parrot.png")

Display the image size.

In [5]:
[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.

In [6]:
plt.figure(figsize=(7, 7))
plt.imshow(imrgb)
plt.show()

It might be useful to convert the color image to gray levels. 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$$

In [7]:
imgray = np.sum(imrgb,2)/3
plt.figure(figsize=(7, 7))
plt.imshow(imgray,cmap='gray',vmin=0,vmax=1)
plt.show()

Now use a custom weighted averaging of the three channels, that reflects better human perception: $$I_{gs}=0.21R+0.72G+0.07B$$

In [8]:
imgray2 = 0.21*imrgb[:,:,0] + 0.72*imrgb[:,:,1] + 0.07*imrgb[:,:,2]
# you can also do: 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()

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.

  1. We compute the histogram of the image imrgb
In [9]:
imhisto, bins = np.histogram(imgray, range=(0,1), bins = 256)
imhisto      = imhisto/np.sum(imhisto)
  1. We now compute the corresponding cumulative histogram thanks to the function np.cumsum()which cumulates the values of a vector from left to right.
In [10]:
imhistocum = np.cumsum(imhisto)
  1. We display the image, histogram and cumulative histogram
In [11]:
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()

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.$$

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.

In [12]:
imeq = imhistocum[np.uint8(imgray*255)]

We can now display the resulting equalized image.

In [14]:
plt.figure(figsize=(5, 5))
plt.imshow(imeq,cmap = 'gray',vmin=0,vmax=1)
plt.show()