Lab on basics in image processing¶
Introduction¶
The objective of this lab is to show how the notions we discovered in the monodimensional case -- that is for signals, can be extended to the two dimensional case. This also enable to have a new look at these notions and perhaps contribute to stenghten their understanding.
In particular, we will look at the problems of representation and filtering, both in the direct (spatial) domain and in the transformed (spatial frequencies) domain. Next we will look at the problems of sampling and filtering.
Within Python, the modules scipy.signal and scipy.ndimage will be useful.
L'objectif de ce laboratoire est de montrer comment les notions que nous avons découvertes dans le cas monodimensionnel - c'est-à -dire pour les signaux, peuvent être étendues au cas bidimensionnel. Cela permet également d'avoir un nouveau regard sur ces notions et peut-être contribuer à renforcer leur compréhension.
En particulier, nous examinerons les problèmes de représentation et de filtrage, à la fois dans le domaine direct (spatial) et dans le domaine transformé (fréquences spatiales). Ensuite, nous examinerons les problèmes d'échantillonnage et de filtrage.
Dans Python, les modules scipy.signal etscipy.ndimage seront utiles.
In order to facilitate your learning and work, your servant has prepared a bunch of useful functions, namely:
rect2 -- returns a two dimensional centered rectangle
bandpass2d -- returns a 2D-bandpass filter (in the frequency domain)
showfft2 -- display of the 2D-Fourier transform, correcly centered and normalized
mesh -- display a ``3D'' representaion of an objet (Ã la Matlab (tm))
To read an image file, you will use the function imread.
To display an image in gray levels, you may use
imshow(S,cmap='gray',origin='upper')
You may either display your graphics inline (it is the default) or using external windows; for that call %matplotlib. To return to the inline mode, use
%matplotlib inline.
Afin de faciliter votre apprentissage et votre travail, votre serviteur a préparé un tas de fonctions utiles, à savoir:
    Rect2 - renvoie un rectangle bidimensionnel centré     Bandpass2d - retourne un filtre passe-bande 2D (dans le domaine de la fréquence)     Showfft2 - affichage de la transformée 2D-Fourier, correctement centré et normalisé     Mesh - affiche une représentation 3D d'un objet (à la Matlab (tm))
Pour lire un fichier image, vous utiliserez la fonction imread.
Pour afficher une image en niveaux de gris, vous pouvez utiliser
Python         Imshow (S, cmap = 'gris', origine = 'supérieur') ......
Vous pouvez afficher vos graphiques en ligne (c'est la valeur par défaut) ou à l'aide de fenêtres externes; Pour cet appel «% matplotlib». Pour revenir au mode Inline, utilisez
% Matplotlib inline.
#This was for Python v. 2.7
#from __future__ import division
#from __future__ import print_function
import matplotlib.pyplot as plt
import numpy as np
from numpy import abs, where, sin, cos, ones, zeros, log, real, imag
from scipy import ndimage as ndi
from scipy import signal as sig
from numpy.fft import fft, ifft, fft2, ifft2, fftshift
from scipy.ndimage import convolve as convolvend
# For 3D representations
from mpl_toolkits.mplot3d.axes3d import Axes3D
# change the dpi for a better vizualisation of images
import matplotlib
savedpi=matplotlib.rcParams['savefig.dpi']
%matplotlib inline
Listing of offered utility functions
# from defs_tp_img import *
#
## Definitions
# ============
def isodd(num):
return num & 1 and True or False
# or num%2==1
##
def rect2(L,N):
""" Returns a (2L+1)x(2L+1) rectangle centered over a
NxN matrix"""
x=np.zeros((N,N))
center=(round((N+1)/2),round((N+1)/2)) if isodd(N) else (round(N/2+1),round(N/2+1))
x[center[0]-L:center[0]+L+1,center[1]-L:center[1]+L+1]=1
return x
##
def bandpass2d(position,half_width,total_dim):
""" Returns the frequency response of a bandpass filter, correctly centered on
`position`, of half-width `half_width` and symmeterized in frequency domain
(zero at the center of the resulting matrix)
Parameters
----------
position: array
center-frequency of the band-pass response
(in samples -- the corresponding frequencies are 2*position/total_dim)
half_width: integer
half frequency band
total_dim: integer
image dimension
Author : jfb"""
N=total_dim
center=(round((N+1)/2),round((N+1)/2)) if isodd(N) else (round(N/2+1),round(N/2+1))
#
H=zeros((N,N))
H[center[0]+position[0]-half_width:center[0]+position[0]+half_width+1,center[1]+position[1]-half_width:center[1]+position[1]+half_width+1]=1
position=-np.array(position) # symmetric frequencies
H[center[0]+position[0]-half_width:center[0]+position[0]+half_width+1,center[1]+position[1]-half_width:center[1]+position[1]+half_width+1]=1
return H
## 2D fequency representation
def showfft2(X,Zero='uncentered', Freq='normalized', num=None,cmap='hot'):
"""Display an image, in the 2D Fourier domain
The representation is centered and in normalized frequencies.
**/!\\ *showfft2* **does not compute the Fourier transform**.
The user is responsible for supplying correct data.
Parameters:
-----------
X : 2D data in the Fourier domain
Num: Figure number(optional)
Zero='uncentered' (default), the input Fourier transform are supposed uncenterd and a
fftshift is applied
Freq='normalized' (default); otherwise axis are in samples
``Author: jfb``
"""
N,M=np.shape(X)
plt.figure(num)
if Zero!='centered':
X=fftshift(X)
if Freq=='normalized':
extent=(-0.5, 0.5, -0.5, 0.5)
else:
extent=(-N/2, N/2, -N/2, N/2)
im=plt.imshow(X,aspect='auto',origin='lower', extent=extent, cmap=cmap)
if Freq=='normalized':
ticks=np.arange(-0.5,0.6,0.1)
plt.xticks(ticks)
plt.yticks(ticks)
plt.colorbar()
def mesh(obj3D,numfig=None,subplt=(1,1,1),cmap='hot'):
"""
Emulates Matlab's mesh
Author: jfb
"""
fig = plt.figure(numfig)
ax = fig.add_subplot(subplt[0],subplt[1],subplt[2], projection='3d')
m,n=np.shape(obj3D)
mm=np.arange(m)
nn=np.arange(n)
X,Y=np.meshgrid(nn,mm)
ax.plot_surface(X, Y, obj3D, rstride=max([round(m/50), 1]),
cstride=max([round(n/50), 1]), cmap=cmap)
def saltpepper(percent=85,shape=None):
s=np.random.random(size=shape)
d=np.where(s>percent/100)
out=np.ones(shape=shape)
out[d]=0
return out
Frequency representation -- Filtering in the frequency domain¶
A pretty sine wave¶
- Generate a sine wave $f(x,y) = \sin(2\pi (f_x x+f_y y))$ , with $f_x$ an $f_y$ between 0.02 Ã 0.2, on $N\times N$ points (e.g. $N=512$). For the implementation, you may use a
forloop, or a double comprehension list, or even the smart functionfromfunction()and an anonymouslambdafunction. Display the result usingimshow.
Note that for a matrix, usually the first index denotes the rows and the second the columns. Therefore modifying something with respect to x results in a variation along the rows of the matrix (usually the "y" axis)..
- Générer une onde sinusoïdale $ f (x, y) = \ sin (2 \ pi (f_x x f_y y)) $, avec $ f_x $ an $ f_y $ entre 0.02 à 0.2, sur $ N \ times N $ points Par exemple $ N = 512 $). Pour la mise en œuvre, vous pouvez utiliser une boucle
for, une double liste de compréhension, ou même la fonction intelligentefromfunction ()et une fonction anonymelambda. Afficher le résultat en utilisantimshow.
Notez que pour une matrice, généralement le premier indice indique les lignes et la seconde les colonnes. Par conséquent, modifier quelque chose par rapport à x entraîne une variation le long des lignes de la matrice (habituellement l'axe \"y\").
?np.fromfunction
fx,fy= 0.02, 0.01
#
# DO IT YOURSELF...
#
# Z=
#
#
#plt.imshow(Z,cmap='hot',origin='upper')
#_=plt.xlabel('y')
#_=plt.ylabel('x')
- Transform your program into a function
pltsinwhich takes the two frequencies as inputs. Then implement a interactive version, using the linesThen play with the parameters, and look at what happens when you change the values of the frequencies.from IPython.html.widgets import interact, interactive, fixed interact(pltsin,fx=[-0.1,0.1,0.01], fy=[-0.1,0.1,0.01])
from ipywidgets import interact, interactive, fixed
def pltsin(fx,fy):
# DO IT YOURSELF...
pass #
#
interact(pltsin,fx=[-0.1,0.1,0.01], fy=[-0.1,0.1,0.01])
- Compute the 2D Fourier transform of $f(x,y)$ (using function fft2) and display the modulus (function abs) of the result, via
showfft2. What are the significations of the axes? What are the spatial frequencies of the sine wave. Experiment with several values of $f_x, f_y$.
# We change the frequencies for ease of visualization in the frequency domain
fx,fy= 0.1, 0.2
#S=
# Display
#
plt.title('2D Fourier transform')
plt.xlabel('Frequencies $f_y$')
plt.ylabel('Frequencies $f_x$')
- Show theoretically, and then check it numerically, that the 2D Fourier transform can be obtained as the succession of two 1-dimensional transforms, applied respectively to the rows and then to the columns (or vice versa). You will take advantage of the fact that the standard
ffthas a parameter axis which is the axis over which thefftis actually computed.
# DO IT!
#
plt.title('2D-Fourier transform as the succession of two 1D transforms')
plt.xlabel('Frequencies $f_y$')
plt.ylabel('Frequencies $f_x$')
Playing with Barbara -- Filtering in the frequency domain¶
Rappelle-toi Barbara
Il pleuvait sans cesse sur Brest ce jour-lÃ
Et tu marchais souriante
Épanouie ravie ruisselante
Sous la pluie (Prévert, 1946)
- Load the image Barbara, using
plt.imread('barbara.png')and display it (in gray levels).
matplotlib.rcParams['savefig.dpi'] = 2*savedpi # double dpi <-- This is to increase resolution
#
#%% reading and displaying the image
#B=
- Display the corresponding frequency representation
showfft2, in log scale (simply takelog(abs())!)
matplotlib.rcParams['savefig.dpi'] = savedpi #restore dpi
# Do it!
plt.title("Barbara's 2D Fourier transform")
- Filter this image using a low-pass filter with a rectangular frequency response (use the function
rect2), for rectangles of half-width 40, 80 ,100. Display the different resulting images, as well as the differences with the initial image. Observations, conclusions.
L=40
#
#
#
#
plt.title('Fourier transform of the filtered image')
The corresponding results in the spatial domain are:
matplotlib.rcParams['savefig.dpi'] = 2*savedpi
#And in the spatial domain
#
#
plt.title("Filtered image (spatial domain)")
# Also experiment with L=80, L=100
# Differences :
#
#
plt.title("Difference between the original image and its low-pass filtered")
- design a frequency response that kills selectively the frequencies around points (46,54) and (-70,79), e.g. on a neighborhood of $\pm$ 10 points. In order to do that, you may use the function
bandpass2d, and you will create a frequency rejector by1-bandpass2d. Look again at the Fourier transform of Barbara, but wit the axes in points, so as to understand what you do. Apply the filter to the initial image and look at the result in the spatial domain. Try to understand the modifications. In particular, look at the tablecloth. Also look at the difference image.
matplotlib.rcParams['savefig.dpi'] = savedpi # restore initial dpi
# Rejector filter
#
#
#
plt.title('Bandpass filter (frequencies)')
matplotlib.rcParams['savefig.dpi'] = savedpi
#
#
#
plt.title("FT of the filtered image")
plt.xlabel('y')
# in the spatial domain we than have
matplotlib.rcParams['savefig.dpi'] = 2*savedpi
#B_filtered
#
#
plt.title("Filtered image by the frequency rejector")
# Differences:
#
#
plt.title("Difference between the initial image and its filtered")
Filtering by convolution¶
The function that will be useful in this part is the function convolve from the module scipy.ndimage (get it by ndi.convolve if ndimage has been imported under the name ndi). We will alo use a np.random.normal(size=(...)) for adding gaussian noise (with typically a scale 0.1); or saltpepper for a salt and pepper noise.
begin by implementing a convolution in two dimensions, by studying the following code:
h=ones((2*ll+1,2*ll+1)) # h the impulse response for m in range(ll,M-ll): for n in range(ll,N-ll): B_filtered[m,n]=sum(sum(h*B[m-ll:m+ll+1,n-ll:n+ll+1]))
Compute a low-pass filtering of Barbara (h constant over an half width of 3 to 10), and examine the result.
(N,M)=np.shape(B)
#
#
#
plt.imshow(B_filtered,cmap='gray',origin='upper')
As we see, it is very simple. Nothing but a sum of products. A better implementation would take care of edge effects.
- Do this filtering again, but using the function
ndi.convolve. Examine the effect of the filtering with an additive gaussian noisenp.random.normal(size=(512,512)with scale 0.1 or a salt and pepper noise. Check that the filtering is indeed a low-pass filter by looking at its transfer function -- use a zero-padding when computing the Fourier transformfft2(h,s=(1000,1000)). You can also use themeshfunction for the representation.
B=plt.imread('barbara.png')
#
#
#
imshow(B_filtered,cmap='gray',origin='upper')
title('Filtered image')
Transfer function
# Lets us look at the transfer function associated with this impulse response..
#
#
#
plt.title('Transfer function of the averaging filter')
On the Barbara image, or on the image of cells, or of bacterias, test a Prewit or a Sobel gradient, with impulse responses
dx=np.array([[1.0, 0.0, -1.0],[1.0, 0.0, -1.0],[1.0, 0.0, -1.0],]) # dx=np.array([[1.0, 0.0, -1.0],[2.0, 0.0, -2.0],[1.0, 0.0, -1.0],]) #Sobel dy=np.transpose(dx)
applied over the two directions (x,y), by ndi.convolve and build the gradient magnitude image.
NB: if Dx and Dy are the gradient images obtained in directions x, y, then the gradient magnitude image is $\sqrt{Dx^2+Dy^2}$.
The ndimage module contains many predefined filters, such as scipy.ndimage.filters.sobel. However, we use here the direct convolution, for pedagogical purposes, instead of these functions.
# Computation of gradient images and magnitude image
#
#
# Sobel Filter
#dx=np.array([[1.0, 0.0, -1.0],[2.0, 0.0, -2.0],[1.0, 0.0, -1.0],])
#dy=np.transpose(dx)
#
#
#
### Prewitt Filter
#dx=np.array([[1.0, 0.0, -1.0],[1.0, 0.0, -1.0],[1.0, 0.0, -1.0],])
#dy=np.transpose(dx)
plt.figure()
#
plt.title('Image filtered using a Sobel gradient over rows')
plt.figure()
#
plt.title('Image filtered using a Sobel gradient over columns')
plt.figure()
#
plt.title("Sobel's Gradient magnitude")
Sobel's transfer function
# Sobel's transfer function
#
#
idem with Barbara
B=plt.imread('barbara.png') # ou cellules.png
#
#
#
#
#
plt.title("Sobel's Gradient magnitude")