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.

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

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

In [3]:
imrgb = plt.imread(path+'parrot.png')


Display the image size:

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

In [6]:
imgray = np.mean(imrgb,2)
plt.figure(figsize=(7, 7))
plt.imshow(imgray,cmap='gray')
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$$

In [7]:
# to do: question 1

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.

a. We compute the histogram of the image imgray:

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

In [9]:
imhistocum = np.cumsum(imhisto) 

c. We display the image, histogram and cumulative histogram:

In [10]:
values = (bins[1:]+bins[:-1])/2
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(15, 5))
axes[0].imshow(imgray,cmap='gray')
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 [11]:
imeq = imhistocum[np.uint8(255*imgray)]

We can now display the resulting equalized image:

In [12]:
plt.figure(figsize=(7, 7))
plt.imshow(imeq,cmap = 'gray')
plt.show()

2. Compute and plot also the corresponding histograms and cumulative histograms of the equalized image.

In [13]:
# to do: question 2

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.

In [14]:
# to do: question 3

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

We start by reading our two images $u$ and $v$:

In [15]:
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$:

In [16]:
index_u = np.argsort(u,axis=None)
u_sort = u.flatten()[index_u] #np.sort(u,axis=None)
index_v = np.argsort(v,axis=None)
v_sort = v.flatten()[index_v] #np.sort(v,axis=None)
uspecifv = np.zeros(nrowu*ncolu)
uspecifv[index_u] = v_sort
uspecifv = uspecifv.reshape(nrowu,ncolu)


We can now display the result.

In [17]:
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(15, 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()