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 = './'
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')
In [4]:
# we can first show imrgb itself, i.e. the NumPy array of values :
print(imrgb)
[[[0.40784314 0.3647059  0.25490198]
  [0.4117647  0.36862746 0.25882354]
  [0.41568628 0.37254903 0.25490198]
  ...
  [0.34117648 0.29803923 0.18039216]
  [0.34117648 0.29803923 0.18039216]
  [0.34117648 0.29803923 0.18039216]]

 [[0.4117647  0.36862746 0.25882354]
  [0.4117647  0.36862746 0.25882354]
  [0.4117647  0.36862746 0.2509804 ]
  ...
  [0.34117648 0.29803923 0.18039216]
  [0.34117648 0.29803923 0.18039216]
  [0.34117648 0.29803923 0.18039216]]

 [[0.4117647  0.36862746 0.2509804 ]
  [0.4117647  0.36862746 0.2509804 ]
  [0.41960785 0.3647059  0.2509804 ]
  ...
  [0.34509805 0.3019608  0.19215687]
  [0.34509805 0.3019608  0.19215687]
  [0.34117648 0.30588236 0.19215687]]

 ...

 [[0.4392157  0.41960785 0.3019608 ]
  [0.4392157  0.41960785 0.3019608 ]
  [0.44313726 0.42352942 0.30588236]
  ...
  [0.30588236 0.28235295 0.19607843]
  [0.30588236 0.28235295 0.1882353 ]
  [0.3019608  0.2784314  0.18431373]]

 [[0.4392157  0.41960785 0.3019608 ]
  [0.44313726 0.42352942 0.30588236]
  [0.44705883 0.42745098 0.30980393]
  ...
  [0.3137255  0.28627452 0.18431373]
  [0.30980393 0.28235295 0.18039216]
  [0.30588236 0.28235295 0.18039216]]

 [[0.44313726 0.42352942 0.30588236]
  [0.44313726 0.42352942 0.30588236]
  [0.44705883 0.42745098 0.30980393]
  ...
  [0.3137255  0.28627452 0.18431373]
  [0.31764707 0.2901961  0.1882353 ]
  [0.30980393 0.28235295 0.17254902]]]


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()
In [7]:
# Let's do the same with a very simple gray level image with only 4 pixels, that we define by hand :
test = np.array([[0,100],[150,200]])

# First we just print the matrix itself :
print(test) 

# Then we use imshow function to display it as an image
plt.figure(figsize=(7, 7))
plt.imshow(test, cmap='gray', vmin=0, vmax=255)
plt.show()
[[  0 100]
 [150 200]]

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

In [9]:
imgray2 = np.average(imrgb, 2, weights=[0.21,0.72,0.07])
# autre solution :
imgray2 = 0.21*imrgb[:,:,0] + 0.72*imrgb[:,:,1] + 0.07*imrgb[:,:,2]
plt.figure(figsize=(7, 7))
plt.imshow(imgray2, cmap='gray', vmin=0, vmax=1)
plt.show()
In [10]:
# The two images imgray and imgray2 look very similar, let's check that they are indeed different :

# First we check that the difference of gray level values are non zero :
print(imgray-imgray2)

# We can also display the difference as an image:
plt.figure(figsize=(7, 7))
plt.imshow(imgray-imgray2, cmap='gray')
plt.show()
[[-0.02359477 -0.0235948  -0.02566013 ... -0.02566013 -0.02566013
  -0.02566013]
 [-0.0235948  -0.0235948  -0.02566016 ... -0.02566013 -0.02566013
  -0.02566013]
 [-0.02566016 -0.02566016 -0.02317649 ... -0.02359477 -0.02359477
  -0.02559477]
 ...
 [-0.02856213 -0.02856213 -0.0285621  ... -0.01981699 -0.02188236
  -0.02188236]
 [-0.02856213 -0.0285621  -0.0285621  ... -0.02346405 -0.02346405
  -0.02394772]
 [-0.0285621  -0.0285621  -0.0285621  ... -0.02346405 -0.02346405
  -0.02552941]]

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 [11]:
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 [12]:
imhistocum = np.cumsum(imhisto) 

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

In [13]:
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()