CS代考 Fundamentals of Image Processing

Fundamentals of Image Processing
Book of Exercises Part 1: Basics
Master 1 Computer Science – IMAge & VCC Sorbonne Universit ́e
Year 2020-2021

Fundamentals of Image Processing – Exercises
IMAge & VCC
Contents
1. Basic Operations
2. Fourier Transform
3. Digitization
4. Image Filtering
5. Filtering and Edge Detection Appendix 1: Dirac delta function Appendix 2: Sampling
Practical works
3 5 7
10 15 19 20 21
Sorbonne Universit ́e
2
Master 1 in Computer Science

Fundamentals of Image Processing – Exercises IMAge & VCC
1. Basic Operations Exercise 1 — Image representations
1. What are the two main methods to represent an image in a computer? Assuming that the first pixel on an image of size N ×M has index 0, how to access the pixel of coordinates (i, j) according to each of these representations? What are the coordinates of the 4-neighbors of pixel (i, j)? How to access these neighbors?
2. What is the size in bits of an image of size 256 × 256 where each pixel value is encoded on one octet (ignore the potential header of the image)? What is the dynamic of this image?
Exercise 2 — Basic matrix operations on images
The following matrices are representations of two images (pixel values in the range [0, 255]):
100 100 100 100 100 100
100 200 200 100 100 100
I= 100 200 200 100 100 100 J=
100 100 100 100 100 100 100 100 100 100 100 100
100 100 100 100 100 100 100 100 200 200 100 100 100 100 200 200 100 100 100 100 100 100 100 100 100 100 100 100 100 100
1. How would you describe intuitively the difference between these two images?
2. Calculate the pixel-wise sum of I and J. What do you observe?
3. Calculate the pixel-wise absolute difference |I−J| of the two images. What do you observe?
4. How is it possible to obtain an image of the same size as I, with only a subpart of I, everything else being black?
Exercise 3 — Histogram
Let us now consider the image represented as:
10 2 1 1 10 2 6 7554423 I=7203218 6320303 5327035
1. What is the minimal grey level kmin in I? the maximal grey level kmax? Derive the dynamics L of I.
2. What is the definition of the histogram? Describe an algorithm to compute the histogram of an image of size N × M encoded in a vector f . Provide the corresponding Python code.
Sorbonne Universit ́e 3 Master 1 in Computer Science

Fundamentals of Image Processing – Exercises IMAge & VCC
3. Draw the histogram of I.
4. Write an algorithm to be applied on the histogram to invert the grey levels of the image.
Can the inverted image be retrieved from the modified histogram? Explain.
5. Write an algorithm to be applied on the image to binarize the image using two values k1 and k2, according to a threshold value S (Ib = k1 if I ≤ S and Ib = k2 otherwise). Draw the histogram of the binarized image.
6. What is a cumulative histogram? How can it be computed and drawn?
7. We want to change the grey level interval [kmin,kmax] into [I1,I2]. Assuming that a linear normalization is performed, write the formula providing the new value of a pixel. Provide the corresponding Python code.
8. Calculate the stretched histogram to the interval [2, 30], and to the interval [2, 8].
9. Histogram equalization consists in applying a transformation providing a flat histogram.
In the continuous case, this is expressed as:
with p′(x) = 1 . kmax −kmin
􏰘 kmax kmin
􏰘 kmax kmin

p(x)dx =
p (x)dx (1)
For the image I above, prove that applying Equation (1) in the discrete case amounts to modify a grey level k into k′ = Int 􏰅 kmax Hc(k)􏰆, with:
N×M
• N the number of rows in the image;
• M the number of columns in the image;
• Hc(k) the value at k of the cumulative histogram; • Int the function returning the closest integer.
Sorbonne Universit ́e 4 Master 1 in Computer Science

Fundamentals of Image Processing – Exercises IMAge & VCC
2. Fourier Transform
Reminder on the continuous Fourier Transform (FT)
Let us recall the definitions of 1D and 2D Fourier Transforms.
1. 1D continuous FT: let x(t), t ∈ R, be a continuous signal. The Fourier transform of x(t),
denoted by X(f), with f ∈ R, is defined􏰘as: +∞
X(f) = x(t)e−2iπftdt (2) −∞
2. 2D continuous FT: let x(t, u) be a continuous signal depending on two variables, t ∈ R, u ∈ R. The Fourier transform of x(t, u), denoted by X(f, g), with f, g ∈ R, is defined as:
􏰘 +∞ 􏰘 +∞
X(f,g) = x(t,u)e−2iπ(ft+gu)dtdu (3)
−∞ −∞
1D continuous FT 2D continuous FT x(t) X(f) x(t,u) X(f,g)
linearity x(t)+λy(t) X(f)+λY(f) x(t,u)+λy(t,u) X(f,g)+λY(f,g)
translation x(t−t0) X(f) e−2iπft0 contraction x(αt) 1 X(f )
Table 1: Main properties of continuous FT.
Exercise 4 — Functional properties of continuous FT
Some useful properties of FT are summarized in Table 1. This exercise aims at proving some of them.
1. Prove that computing the 2D FT is equivalent to compute successively two 1D FT, i.e. X(f, g) = F T [Z(f, u)]. Explain the function Z(f, u) and provide its mathematical expres- sion.
2. Prove that, if a 2D function x(t, u) is separable (can be written as the product of two 1D functions as x(t, u) = z(t) × k(u)), then we have: X(f, g) = Z(f) × K(g), with Z(f) = F T [z(t)] and K(g) = F T [k(u)].
3. Prove the links between convolution (denoted by ⋆) and product: FT[x⋆y(t)]=X(f)×Y(f)andFT[x×y(t)]=X⋆Y(f). Recallthattheconvolutionof two functions x(t) and y(t) is defined as:
􏰘 +∞ −∞
|α| α
convolution x⋆y(t) X(f)×Y(f) x⋆y(t,u) X(f,g)×Y(f,g)
x(t−t0,u−u0) x(αt,βu)
X(f,g) e−2iπ(ft0+gu0) 1 X(f , g )
|α||β| α β product x(t)×y(t) X ⋆Y(f) x(t,u)×y(t,u) X ⋆Y(f,g)
z(t) = x ⋆ y(t) = Sorbonne Universit ́e 5
x(τ)y(t − τ)dτ (4) Master 1 in Computer Science

Fundamentals of Image Processing – Exercises IMAge & VCC 4. Prove the rotation property of the 2D FT:
F T [x (t cos θ + u sin θ, −t sin θ + u cos θ)] = X (f cos θ + g sin θ, −f sin θ + g cos θ) Exercise 5 — Computing the FT of some usual signals
1. From Equation (2), compute the FT of the function x(t) = Rect 􏰁 Tt 􏰂, where T is a constant and Rect(t) is defined as:
􏰋 1 if |t| ≤ 1
Rect(t) = 2 (5)
Derive the FT of the corresponding 2D function x(t, u) = Rect 􏰁 Tt 􏰂 × Rect 􏰁 Tu 􏰂. 2. Let δ(t) denote the Dirac delta function (see Appendix 1):
δ(t)=􏰋 0 if t̸=0 ∞ otherwise
Let us recall that FT [δ(t − t0)] = e−2iπft0 and FT 􏰃e2iπf0t􏰄 = δ(f − f0). Derive:
(a) the FT of x(t) = cos (2πf0t) and x(t) = sin (2πf0t); (b) the FT of a 2D Dirac delta function δ(t − t0, u − u0); (c) theFTofsθ(x,y)=Acos[2πfo(xcos(θ)+ysin(θ))].
0 otherwise
Sorbonne Universit ́e 6 Master 1 in Computer Science

Fundamentals of Image Processing – Exercises IMAge & VCC
3. Digitization Reminder on digital signals and images
Signals and images that are processed and analyzed computationally are digital. This impacts their frequential analysis in three main ways:
1. The analysis is performed on a limited temporal or spatial domain: how are the FT of a continuous signal with infinite band and the FT of a band-limited continuous signal related? (see Exercise 7).
2. Sampling of a continuous signal, i.e. choice of a countable set of values: what is the FT of the sampled signal? Is some information lost? in which cases? (see Exercise 6).
3. Quantization of the signal (see lecture) & Discrete Fourier Transform (DFT) (see Exer- cise 8) .
Exercise 6 — Sampling a continuous signal
The ideal sampling consists in selecting in the continuous signal a set of regularly spaced values, forming a discrete (or digital) signal xe(t).
This writes (see Appendix 2):
􏰒+∞
xe(t) = x(t)δ(t − kTe) = x(t) Te (t) (6)
k=−∞
where Te is the sampling period, and Te (t) is called Dirac comb (or impulse train), built from
Dirac delta functions for the given period Te (see Appendix 2).
1. Prove that the FT of the sampled signal xe(t) can be written as:
1􏰒+∞ k
Xe(f) = T X(f − T ) (7)
e k=−∞ e
Indication: Use the properties listed in Appendix 2.
Interpretation: The FT of the sampled signal xe(t) corresponds to the FT of the contin-
uous signal x(t) periodised with period fe = 1 . Te
2. Reconstruction of the continuous signal from the digital signal.
(a) Propose a method in the frequential domain to recover the FT X(f) of the continuous signal from the FT Xe (f ) of the sampled signal. Which condition on X (f ) and fe should hold for this method be possible? This condition is known as Shannon theorem.
(b) Under the Shannon condition, it is possible to reconstruct a continuous signal x(t) from the samples xe(kTe), by inverting Equation (7). Prove that the reconstruction writes:
x(kTe)sinc(πfe(t − kTe)) (8) Master 1 in Computer Science
+∞ +∞
xr(t) = 􏰒 x(kTe)sin(πfe(t − kTe)) =
k=−∞ πfe(t − kTe) k=−∞
Sorbonne Universit ́e 7
􏰒

Fundamentals of Image Processing – Exercises IMAge & VCC
3. Spectral overlapping.
If the Shannon condition is not satisfied (because fe is too low or because it is impossible to satisfy it), what is the consequence on the spectrum of the sampled signal?
Example: sinusoidal signal
Let c(t) = A cos(2πf0t) be a sinusoidal signal.
(a) Let this signal be sampled at frequency fe, yielding ce(k), k ∈ {1;N}. Is there a range of sampling steps that allow for the reconstruction of the signal? In that case, give this range.
(b) The signal is now sampled at frequencies 4f0 and 2f0. Draw the spectrum of the sampled signals at these frequencies.
(c) Now c(t) is sampled at frequency fe = 32f0. Draw the spectrum of ce(t). What is observed? What is the reconstructed signal if Equation (8) is applied?
Exercise 7 — Signals with limited support.
Let x(t) be a continuous signal, and x (t) the signal which is equal to x(t) in the interval L
􏰃􏰄􏰋􏰁􏰂
−L2 ; L2 , L ∈ R, and 0 elsewhere, i.e. xL(t) = x(t)Rect Lt , where Rect(t) is again:
1 if |t| ≤ 1
Rect(t) = 2 (9)
0 otherwise
1. Write XL(f) = FT [xL(t)] as a function of X(f).
2. Compute XL(f) for x(t) = Acos(2πfot). Draw |XL(f)|.
3. Is it theoretically possible to recover X(f) from XL(f)? For the example, propose a method that guarantees that XL(f) will suffer only a low degradation with respect to X(f). Indication: The amplitude of secondary lobes of the sinus cardinal function is low in comparison to the amplitude of the main lobe. It is therefore possible to consider that the amplitude of sinc(t) becomes 0 after k secondary lobes.
4. Bonus: Same question in 2D: compare x(t, u) and xL (t, u) = x(t, u)Rect 􏰁 Lt 􏰂 Rect 􏰁 Lu 􏰂, for x(t, u) = A cos [2πfo (t cos(θ) + u sin(θ))].
Exercise 8 — DFT and frequency resolution
The aim of this exercise is to study the three discretization steps on a particular example. Let c(t) = Acos(2πf0t)Rect􏰁Lt 􏰂 (see Exercise 7). Here we set L = 2T0. 􏰑
1. Let sample c(t) at frequency fe = 4f0. Draw the sampled signal ce(t) = k c(t)δ(t − kTe), and the digital signal c(k).
2. Draw the FT of the sampled signal.
Sorbonne Universit ́e 8 Master 1 in Computer Science

Fundamentals of Image Processing – Exercises IMAge & VCC
3. Let us recall that the DFT consists in sampling the FT of the sampled signal in the
interval [−fe ; fe [. Draw the DFT X(h) of c(k). Explain why the secondary lobes of the 22
sinus cardinal do not appear.
4. Let us now sample c(t) at frequency fe = 3.7f0. Compute X(h) and provide a graphical representation. What do you conclude?
5. Zero-Padding: let p(k) be the following discrete signal with 2Nvalues:
p(k) = 􏰋 x(k) if k ∈ {0; N − 1} (10)
0 if k ∈ {N ; 2N − 1}
Compute the DFT P(h) = DFT [p(k)], and compare it with X(h). What do you conclude
on the effect of zero-padding?
Sorbonne Universit ́e 9 Master 1 in Computer Science

Fundamentals of Image Processing – Exercises IMAge & VCC
4. Image Filtering Reminder on spatial and frequential filtering
Image filtering consists in processing the frequential content of an image in order, for example, to:
1. keep only the low frequencies (low-pass filter), e.g. to attenuate the noise;
2. keep only the high frequencies (high-pass filter), e.g. to detect contours (see exercise session 5);
3. select a given range of frequencies (band-pass filter).
Let x(n, m) denote the image value at pixel (n, m). Frequential filtering consists in multiplying
the spectrum (i.e. the FT) X(f, g) by a function H(f, g), yielding the filtered spectrum XF (f, g):
XF (f, g) = X(f, g)H(f, g) (11)
This operation in the frequency domain corresponds to a convolution operation (⋆) in the spatial domain of x(n, m) by the inverse FT h(n, m) of H(f, g):
xF (n, m) = x ⋆ h(n, m) (12)
The function h(n, m) is called impulse response of the filter, and corresponds to the filtered image of an image reduced to a Dirac impulse at (0, 0).
In the digital (discrete) case, if the support of h is finite (of size d, supposed to be odd), the filter is called Finite Impulse Response (FIR) filter, and the filtering operation writes:
d−1
􏰒2 􏰒2
xF (n, m) =
i=− d−1
x(i, j)h(n − i, m − j) (13) j=− d−1
d−1
22
Filtering is then performed according to the following steps:
1. Apply a rotation of angle π of h with respect to its center: h(n, m) ⇒ h(−n, −m) = g(n, m);
2. Center the filter at pixel p;
3. Computer the weighted sum of the image pixels and the filter coefficients g(n, m); 4. This weighted sum is the value at p of the filtered image.
Exercise 9 — Filters
1. Which filter has for impulse response h1(u,v) a 2D Rect (“gate”) function, centered and of size 3?
1 0 -1
2. Writethefilterh2(u,v)= 1 0 -1 asasumofshiftedDiracfunctions.
1 0 -1
Sorbonne Universit ́e
10 Master 1 in Computer Science

Fundamentals of Image Processing – Exercises IMAge & VCC
3. Show that the filter h2 is separable, and express it as h2 = (f1)T f2, where f1 and f2 are 1D convolution kernels.
4. Boundary effect. An image has a bounded support, hence pixels on the boundary do not have neighbors outside the support, which makes the usual application of the convolution formula impossible. What would you propose to compute a filtered value at those pixels?
123
5. Consider the image I = 4 5 6 . Compute I ⋆h2 and (I ⋆f1T)⋆f2. What do you
789
conclude?
6. Generalization: let us consider a separable linear filter h(n, m) = f1 (n) · f2 (m). Prove that the 2D filter can be decomposed into two successive 1D filters.
Exercise 10 — Some filters
Let consider the following 8 × 8 image, with values (grey levels) ranging from 0 to 7:
00000000 05555550 05777750 05777750 05777750 05777750 05555550 00000000
1. Draw the grey value profile of line 3 of this image (lines are numbered from 0 to 7).
111
2. Filter this image using a linear convolution by the filter 19 × 1 1 1 (compute the grey
111
levels of all pixels of the filtered image). Draw the grey value profile of line 3 of this filtered
image. What do you observe?
0 -1 0
3. Same question with the filter -1 4 -1 0 -1 0
4. One wants to filter an image with a Gaussian filter of standard deviation σ = 0.5. Compute the coefficients of the filter. Draw the profile of the centerline of the filter. Same questions with σ = 1.
5. Let us consider that the Gaussian function becomes 0 after 3σ. What is the size of the
Gaussian mask in the spatial domain to be used in order to obtain a filter with cut-off
frequency (in the frequential domain) equal to fe ? 2
Sorbonne Universit ́e 11
Master 1 in Computer Science

Fundamentals of Image Processing – Exercises IMAge & VCC
Exercise 11 — Convolution
Let us apply the convolution operator on a small region of size 20 × 11 (in red in the figure) of the following image. This region contains a part of the spinnaker pole, and just below a rope.
The following matrix I contains the grey levels in this region:
196 196 199 203 200 188 195 201 204 198 200 198 212 188 206 199 202 204 201 204 187 213 199 187 206 188 209 181 213 208 208 194 186 198 218 202 198 194 210 195 195 208 163 128 184 194 187 206 207 223 202 181 149 125 81 207 195 206 217 205 195 174 151 108 70 74 205 198 187 193 190 161 133 94 80 62 210 179 196 205 186 155 119 75 66 139 223 209 204 196 185 115 91 77 60 189 202 208 188
I=208 159 105 95 68 90 222 214 206 199 191 146 100 87 55 151 210 210 200 192 207 201 84 69 64 207 212 212 194 211 195 196 147 60 108 219 194 198 187 188 204 191 214 208 174 218 209 195 203 190 205 192 201 201 192 200 215 201 187 204 195 203 192 210 198 148 207 194 192 184 198 197 188 205 164 161 208 176 191 200 195 191 197 207 180 176 217 198 197 190 201 200 208 195 154 171 198 198 192 193 189 205 198 216 141 192 203 194 199 196 196 217 207 186 146 207 211 188 203 188 180
1. Which are in this image the pixels belonging to the spinnaker pole and to the rope?
2. Low-pass filtering. Let us consider a binomial low-pass filter of size 3 × 3, whose convo- lution mask is:
121 h1=1 2 4 2
16
Sorbonne Universit ́e 12
121
Master 1 in Computer Science

Fundamentals of Image Processing – Exercises IMAge & VCC
Applying this mask on the considered 20 × 11 region provides the following result (not complete):
Complete the result by computing the missing values (in grey), rounded to their integer part.
3. Contrast enhancement. In order to enhance the contrast on the edges of the objects, the following convolution mask is applied:
0 -1 0 h2=1 -1 5 -1
16
0 -1 0 The resulting filtered image (incomplete) is:
Sorbonne Universit ́e 13 Master 1 in Computer Science

Fundamentals of Image Processing – Exercises IMAge & VCC
Complete the result by computing the missing values (in grey).
Sorbonne Universit ́e 14 Master 1 in Computer Science

Fundamentals of Image Processing – Exercises
IMAge & VCC
5. Filtering and Edge Detection
Exercise 12 —
Let I be defined as:
 9 8 9 9 9 9 8 9 8 7   9 8 9 8 8 9 9 8 6 4   8 8 9 8 8 8 8 6 4 3   9 9 9 8 9 8 6 4 3 3 
I=8 9 9 8 8 6 4 3 2 2  9 8 8 8 6 5 3 3 1 1   8 8 8 5 5 3 2 2 1 0   8 6 5 4 2 3 2 1 0 0 
7542221000 4432111000
and corresponding to the following image:
to which Sobel filters will be applied.
Sobel edge detection filter
Sorbonne Universit ́e 15 Master 1 in Computer Science

Fundamentals of Image Processing – Exercises
The impulse responses of Sobel filters (convolution masks) are defined as:
10−1 1 2 1 Sx= 2 0 −2 andSy= 0 0 0
1 0 −1 −1 −2 −1
1. Complete the computation of the horizontal filtering Ix = I ⋆ Sx: 000000000000  0 −4 0 3 −1 1 −2 −1 −3 −10 −5 0   0 −3 1 1 ? 2 1 −4 −10 −13 ? 0  0 −1 2 −1 −3 1 −2 ?−14−11−40  0 1 2 −3 −2 −2 −10 ? −12 −6 −1 0 
I =0 1 1 −3 −4 −7 ? −12 −9 −5 0 0
000000000000 N.B: Padding with copy on the image I is used here.
IMAge & VCC
x0−1−1 −4−8?−13−8−7
 0 −3 ? −8 −11 −8 −9 −6 −6
0−6 ?−10−11−4 −4 −7 −6
 0 −6 ? −10 −9 −2 −2 −7 −5 −1 0 0   0 ? −6 −9 −8 −3 −1 −5 −4 0 00
−7 −1 0  −7 −2 0  −4 −10
2. Complete the computation of the vertical filtering Iy = I ⋆ Sy : 000000000000 0 0 0 −1 −3 −3 0 1 −3 −8 −11 0 0 −3 −1 −1 ? −4 −3 −4 −10 −15 ? 0 0 1 2 1 1 1 −4 ? −14 −11 −6 0 0 1 2 1 0 −2 −8 ? −12 −8 −5 0
I=0 −1 −3 −3 −4 −9 ? −10 −7 y 0 −1 −3 −6 −10 ? −11 −8 −5 0 −5 ? −12 −15 −14 −9 −6 −6 0 −6 ? −14 −13 −10 −6 −5 −6
−7 −8 0 −5 −7 0 −5 −4 0 −4 −1 0
0 −14 ? −8 −7 −6 −6 −5 −3 −1 0 0 0 ? −6 −3 −2 −3 −3 −1 0 0 0 0 000000000000
Sorbonne Universit ́e 16 Master 1 in Computer Science

Fundamentals of Image Processing – Exercises
IMAge & VCC
→−
3. Complete the computation of the norm of the gradient ||G|| =
closest integer):
􏰛
Ix2 + Iy2 (rounded to the
000000000000  0 4 0 3 3 3 2 1 4 13 12 0   0 4 1 1 ? 4 3 6 14 20 ? 0   0 1 3 1 3 1 4 ? 20 16 7 0   0 1 3 3 2 3 13 ? 17 10 5 0 
→− 0134611?1611980 ||G||=0 1 3 7 13 ? 17 11 9 9 7 0  0 6 ? 14 19 16 13 8 8 9 4 0   0 8 ? 17 17 11 7 9 8 6 1 0   0 15 ? 13 11 6 6 9 6 1 0 0  0?8984354000
000000000000
4. Propose a threshold value to binarize the image and extract the contours from ||G||.
Exercise 13 — Suppression of non-maxima
The problem of edges detection can be formulated as the selection of local maxima of the
gradient amplitude in the direction of the gradient. A discrete filter approximating the gradient
is first applied, and then a post-processing step allows suppressing the non-maxima points.
Let G(P ) be the image gradient of a pixel P . One has to check whether P is a local maximum −−−→
(a) (b) Two algorithms are considered to do this:
1. Discretization of orientations. Propose an algorithm to determine the discrete direction to be assigned to point P (considering 8-neighborhood).
2. Sub-pixel interpolation. Propose a method based on the computation of the gradient of the two neighbors that are in the direction of the gradient.
Indication: Determine the coordinates of points Q and R in Figure (b) above, and compute the gradient module using linear interpolation.
Sorbonne Universit ́e 17 Master 1 in Computer Science
−−−→
of ||G(P)|| in the direction of the gradient (Figure (a) below).
→−

Fundamentals of Image Processing – Exercises IMAge & VCC
Exercise 14 — Multi-resolution pyramid and aliasing
We assume a multi-resolution pyramid built from a Gaussian filtering. The aim here is to estimate the parameters of the filter so as to not induce aliasing effect during sub-sampling. It is assumed that the original image (with the best resolution) has been correctly sampled.
Recall that the FT of a Gaussian function is a Gaussian function with inverted standard
deviation (σf = 1 ): σsπ
􏰙 −b2(t2+u2)􏰚 π −π2(f2+g2) TFe =eb2
b2
1. By using a Gaussian filter, is it theoretically possible to cancel totally the aliasing effect?
Explain.
2. Assuming that in practice the Gaussian filter has a cut-off frequency at 3σ, determine the standard deviation of the Gaussian mask to be applied in the spatial domain, so as to guarantee that the convolution I ⋆ G does not introduce aliasing when moving to the next resolution level.
Sorbonne Universit ́e 18 Master 1 in Computer Science

Fundamentals of Image Processing – Exercises
IMAge & VCC
Appendix 1: Dirac delta function
Intuitively, the Dirac delta function can be interpreted as follows: δ(t) = lim 􏰃T1 Rect􏰁Tt 􏰂􏰄. T→0
Exercise 15 — Properties of the Dirac delta function
1. 􏰗+∞δ(t)dt=1. −∞
2. x(t)δ(t − t0) = x(t0)δ(t − t0)
3. x⋆δ(t−t0)=x(t−t0)
4. scaling property: |α|δ(αt) = δ(t)
5. FT [δ(t − t0)] = e−2iπft0 and FT 􏰃e2iπf0t􏰄 = δ(f − f0)
Sorbonne Universit ́e 19 Master 1 in Computer Science

Fundamentals of Image Processing – Exercises
IMAge & VCC
Appendix 2: Sampling
The Dirac distribution δ(t) is defined as: δ(t)=􏰋 0
if t̸=0 otherwise
(14)
Dirac Comb
The Dirac comb is defined as:
+􏰑∞ k=−∞
Sampling model
Te (t) =
δ(t − kTe).
The ideal sampling of a 1D signal x(t) with periode Te is:
Signal FT Signal FT 1 δ(f) δ(t) 1
􏰒+∞ k=−∞

x(t)δ(t − kTe) = x(t)
Exercise 16 — FT of Dirac distribution and other usual functions
xe(t) =
With this model, we have lim xe(t) = x(t) and lim 􏰗 +∞ xe(t)dt = 􏰗 +∞ x(t)dt.
Te (t) Te→0 Te→0 −∞ −∞
(15)
e−2iπf0t
+∞ +∞
δ(f − f0) δ(t − t0) e−2iπft0
􏰑􏰑
T ( t ) = δ ( t − k T ) T1 1 ( f ) = T1 δ ( f − Tk )
k=−∞ T k=−∞ 􏰅􏰆
Rect􏰁 t 􏰂 Tsinc(πfT) f0sinc(πf0t) Rect f T f0
Sorbonne Universit ́e 20 Master 1 in Computer Science

Fundamentals of Image Processing – Exercises IMAge & VCC
Practical works
Practical work materials (Jupyter notebooks and data) are to be retrieved from the course’s web site: https://frama.link/mR68kcMB.
Sorbonne Universit ́e 21 Master 1 in Computer Science

Practical work 1: introduction and image enhancement
• Quick start for Python (10 minutes!) : https://www.stavros.io/tutorials/python/
• Quick start for Numpy : https://numpy.org/devdocs/user/quickstart.html#
• For Matlab users: Numpy is very similar but with some important difference, see
http://mathesaurus.sourceforge.net/matlab-numpy.html.
• Keep in mind that in Python, exception of variable of scalar type, all is reference and affec-
tation is not a copy.
Short introduction to image processing with Python
Help: use the function help() to get information on a Python objet.
Images are stored as arrays that is the default type of the numpy module. Defaut type of array elements is float64 according to the IEEE754 norm. Special float values are defined: infinity (inf) and undefined (nan, not a number), and some numerical constants, such as π.
[ ]:
# import numpy
import numpy as np
# predefined constants
print(np.inf,np.nan,np.pi)
# some values
print( 1., 1e10, -1.2e-3)
Creating an array: several ways.
1. From a list of values (formally any Python iterable object). Elements of an array have the same type, determined by Numpy:
[ ]:
[ ]:
3. From a sequence, prefer arange() from numpy to range() from python. [ ]:
1
V = np.array([1,2,3])
M = np.array([[1,2,3],[4,5,6.]])
print (“V is of type”,V.dtype)
print (“M is of type”,I.dtype)
2. Without values: Numpy has constructors such as empty(), zeros(), ones(). . . Shape should be given (see below). Important: empty() does not initialize array elements.
I = np.zeros((3,4))
print(I)
J = np.empty((4,3))
print(J)
print(np.arange(10))
print(np.arange(0,10,2))
print(np.arange(9,-1,-.5))

Shape of an array
Shape decribes the number of elements for each dimension. A vector is of dimension 1, a matrix is of dimension 2. Superior dimensions are possible. Shape is not size that is the number of elements of an array. Type of shape is always a tuple of integers. With previous example:
[ ]:
[ ]:
[ ]:
[ ]:
print(I.shape, I.size)
print(J.shape, J.size)
print(V.shape, V.size)
An important function/method is reshape() to change the shape of an array. Typical usage of reshape() is to transform a vector into a matrix or reciprocally.
K = np.arange(12)).reshape((3,4))
print(K)
print(np.reshape(K,(12)))
print(K.reshape((2,2,3)))
Elements of an array
Access element by indices: two syntaxe are possible, the first given in the example is prefered. Negative index is possible with the same meanning of Python list.
I = np.arange(12).reshape((3,4))
print(I[1,2])
print(I[0][0])
print(I[-1,0])
Access by group of indices using the operator : allows to extract subarray. General syntaxe is start:end:step and it is very powerfull:
print(‘extract the first line’)
print(I[0,:])
print(I[0,0:])
print(I[0,::])
print(I[0,::1])
print(‘extract center of the array’)
print(I[1:3,1:3])
print(‘extract elements with even indices’)
print(I[::2,::2])
print(‘print the horizontal mirror of an array’)
print(I[:,::-1])
2

Array arithmetic
Operators and functions can be applied to arrays. Mostly, operations are element-wise (i.e. applied element by element). The consequence is arrays should have the same shape. One operand can also be scalar in most of time.
[ ]:
[ ]:
[ ]:
A = np.arange(12).reshape((3,4)) B=2*A+1
C=A+B
D = np.cos(2*np.pi*A/12)
print (D)
print (D**2)
print (D>0)
Array may be viewed as matrix, we can make some linear algebraic manipulation. For example, np.matmul() is the matrix multiplication. It can be used to build matrix from vector. An example, using the transpose operator T.
L = np.arange(1,6).reshape((1,5))
# transpose of L. Warning: C remains a reference to L C = L.T
# This could be better if your want to touch L
C = L.T.copy()
print(“A 5*5 matrix:”)
print(np.matmul(C,L))
print(“A dot product, but result is a matrix:”)
print(np.matmul(L,C))
print(np.matmul(L,C)[0,0])
print(“dot() is prefered with vectors:”)
V = np.arange(1,6)
print(V.dot(V))
print(np.dot(V,V))
Images
We make use of PIL module (https://pillow.readthedocs.io/en/stable/reference/Image.html) to load and write an image and easily converted to Numpy array. Be careful: array type depends on image.
from PIL import Image
# reading an image and convert to array
myimage = np.array(Image.open(‘image.png’))
3

# write an image (alternative format) from an array
Image.fromarray(myimage).save(‘image.jpg’)
Array can be displayed as an image using Matplotlib module. Here a short example:
[ ]:
import matplotlib.pyplot as plt
# minimal example:
plt.imshow(myimage)
plt.show()
# with more controls:
w,h=400,400
plt.figure(figsize=(w/80,h/80)) # optional, to control the size of figure (unit:
􏰀→ pixel)
plt.gray() # optional call to display image using a gray colormap plt.title(‘This is an image’) # optional: add a title plt.axis(‘off’) # optional: remove axes
plt.imshow(myimage)
plt.show()
See also: – https://matplotlib.org/3.1.1/tutorials/introductory/images.html – https://matplotlib.org/gallery/images_contours_and_fields/image_demo.html#sphx-glr- gallery-images-contours-and-fields-image-demo-py).
Exercice 1
In this exercice, we work with image img/moon.tif. If possible give two solutions : one with loops (for, while, . . . ) and one without loops.
1. Write and test a function openImage() getting an image filename as argument and returning the array of pixel values.
[ ]:
[ ]:
from PIL import Image
import numpy as np
def openImage(fname):
“”” str -> Array
(notation above means the function gets a string argument and returns an␣
􏰀→Array object) “””
2. Write and test a function countPixels() getting an array and an integer k as arguments and returning the number of pixels having the value k.
def countPixels(I,k):
“”” Array*int -> int”””
4

3. Write and test a function replacePixels() getting an array and two intergers and replacing pixels having k1value to k2 value and returning the new array. Be aware to not modify I.
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
def replacePixels(I,k1,k2):
“”” Array*int*int -> Array “””
4. WriteandtestafunctionnormalizeImage()gettinganarrayandtwointegersk1andk2and returning an array with elements normalized to the interval [k1, k2].
def normalizeImage(I,k1,k2):
“”” Array*int*int -> Array “””
5. Write and test a function inverteImage() getting an array and returning and arry having inverted pixel values (i.e. the transform k 􏰖→ 255 − k
def inverteImage(I):
“”” Array -> Array “””
6. Write and test a function computeHistogram() getting an array and returning its histogram. Type of histogram can be an array or a list. It is forbidden to use an histogram method from a Python module. Is it possible to compute the histogram without explicitely visiting array pixels?
def computeHistogram(I):
“”” Array -> list[int] “””
# use comments to answer to a verbal question
7. Write and test a function thresholdImage() getting an array I and an integer s and return- ing an array having elements set to 0 if corresponding element of I is lower than s or 255 otherwise.
def thresholdImage(I,s):
“”” Array*int -> Array “””
8. Using previous functions, give a series of instructions to read then to display an image, plot the histogram (one can use plot() or bar() from matplotlib.pyplot module), inverse the image and display it, plot its histogram.
import matplotlib.pyplot as plt ## your code start below
9. Give a series of instructions to read and display an image, plot the histogram, normalize the image to the interval [10, 50], compute the new histogram, display the image and the histogram. Remark: imshow() normalizes image. To avoid this and see the effect of the normalization, use imshow() with parameters vmin=0,vmax=255. Comment the results.
5

10. Same question than 9. remplacing the normalization by a thresholding with parameter s = 127.
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
Exercice 2 – generate images
1. Create the array I 4 by 4 corresponding to the following image:
Black pixels have value 0, white pixels value 255, and grey pixels value 127. Display the image using imshow() and plot the histogram.
2. We want to generate a matrix having random values. Functions rand() and randn() from numpy.matlib module generate array of given shape with random values following respec- tively a uniform distribution on [0, 1[ and a normal distribution. Create an array of shape 512 by 512 having integer elements following an uniform distribution in the set {0, 1, · · · , 255} . We also want to create an array following a gaussian distribution with a mean of 128 and a standard deviation of 16 and with integer values. Display the images and their his- togramms. Discuss the results.
import numpy.matlib
Exercice 3: image manipulation
In this exercice, we work with image img/pout.tif. 1. Read and display this image
2. Examine the histogram. Determine the extrema of the image. What can you say about the quality of this image?
3. Using functions from Exercice 1, write the function histogramEqualization() getting one image, its histogram, applying an histogram equalization and returning the new image. Test this function on pout.tif and discuss the result.
def histogramEqualization(I,h):
“”” Array * (list[int] -> Array “””
6

Practical work 2 : Fourier transform
This practial work is dedicated crete Fourier transform applied
to the study of the dis-
on
the two following
images:
1

and analyze the properties of their spectrum. To this end, we make use of the following functions
provided by the module numpy.fft:
• fft2() to compute the Fourier transform on an image
• fftshift() to center the low frequencies
• abs() (from numpy) to compute the module of a complexe array
In most of cases, high frequencies have lower energy compare to low frequencies. We will use a logarithmic scale by applying log(1 + abs(TF)) to display the spectrum.
[ ]:
import numpy as np
from numpy.fft import fft2,fftshift
from PIL import Image
son = np.array(Image.open(‘img/son.png’))
sonrot = np.array(Image.open(‘img/sonrot.png’))
2

Exercice: properties of Fourier transform applied on natural images
1. Write the following functions:
• computeFT(I) returning the Fourier transform of image I,
• toVisualizeFT(If) returning the centered module of a complex array If (the Fourier trans-
form of an image),
• toVisualizeLogFT(If) similar to the previous function but use a logarithmic scale.
[ ]:
[ ]:
[ ]:
[ ]:
3. Interpretation: discuss the results obtained on thresholded FT module. What property of the Fourier transform is shown ?
4. Write the function blend() getting two images, one float parameter α ∈ [0, 1], calculating αI1 + (1 − α)I2 and returning the result.
def computeFT(I):
“”” Array -> Array[complex] “””
def toVisualizeFT(If):
“”” Array[complex] -> Array[float] “””
def toVisualizeLogFT(If):
“”” Array[complex] -> Array[float] “””
2. Write a series of instructions that
• compute the Fourier transform of son and sonrot,
• compute and display the module using a logarithmic scale,
• threshold the module with a parameter of 1.105 (use the function of TME1)
• display the thresholded spectrum
import matplotlib.pyplot as plt # your code below
def blend(I1,I2,alpha):
“”” Array**2*float -> Array “””
5. Apply the previous function on images son and sonrot and α = 21 , compute the Fourier transform, threshold the module and visualize the result.
6. Compare the latter result with those of question 2: what property of the Fourier transform is shown? What is the behaviour of α in the resulting spectrum?
7. We want to determine the text orientation in image sonrot and produce a new image with horizontal text. Write the function rectifyOrientation() that:
• computes the FT module of image given in parameter and threshold it at 3 × 105,
• from thresholded module determines the main orientation using the function
mainOrientation()
3

• produces the rectified image applying a rotation with a suitable angle using rotateImage() [ ]:
[ ]:
def mainOrientation(I):
“”” Array -> tuple[Iori,float]
return image of orientation (32 bins) and the main orientation (degree)␣ 􏰀→from a Fourier transform module
“””
n, m = I.shape
size = 32
x = np.array(range(size))
ori = np.vstack((np.cos(np.pi*x/size), np.sin(np.pi*x/size))).T
Iori = np.zeros((n, m))
orients = np.zeros((size))
for i in range(1,n+1):
for j in range(1,m+1):
if I[i-1, j-1] > 0:
v = np.array([j-m/2, -i + n/2])
if i > n/2:
v = -v
prod = np.matmul(ori, v)
maxi = prod.max()
if maxi > 0:
imax = np.nonzero(prod == maxi)
Iori[i-1, j-1] = imax[0]
orients[imax] += 1
maxori = np.nonzero(orients == orients.max())[0][0]
return (Iori, 180*maxori/size – 90)
def rotateImage(I,a):
“”” Array*float -> Array
return a rotation of angle a (degree) of image I
“””
return np.array(Image.fromarray(I).rotate(a, expand=True, fillcolor=127)) ####### your code below
8. Experiment rectifyOrientation() on sinrot, and on a rotation of img/port.jpg (using rotateImage()) with various rotation angles.
4

Practical work 3: 2D sampling and aliasing
Properties studied in 1D apply in 2D. The following results can be admitted : – given a regular grid, a sampling of a continuous 2D signal can be modelled as follow:
+∞ +∞
xe(t,u) = x(t,u)e(t,u) with e(t,u) = ∑ ∑ δ(t−kTe,u−lTe)
k=−∞ l=−∞
e is the analog of Dirac comb (also called impulse train) in 2D – spectrum of xe writes:
1 +∞ +∞
Xe(f,g)=Te2 ∑ ∑X(f−kfe,g−lfe)
k=−∞ l=−∞
2D sampling then implies a periodisation of the spectrum for the two dimensions – it is possible to reconstruct the original signal from the sampled signal if 2D Shannon condition is verified (band limited signal) with:
+∞ +∞
xr(t, u) = ∑ ∑ xe(kTe, lTe) sinc(π fe(t − kTe)) sinc(π fe(u − lTe)) (1)
k=−∞ l=∞
so called Shannon interpolation.
Exercice 1: aliasing and windowing of 2D signals
Given the following signal:
sθ(t, u) = A cos(2π f0(t cos θ + u sin θ)) Here an example with θ = π4 :
The goal of this exercice is to study the limit conditions of sampling of this image in order to avoid aliasing.
[ ]:
import numpy as np
from numpy.fft import fft2, fftshift
import matplotlib.pyplot as plt
# for interactive ploting, see surf() below %matplotlib notebook
from matplotlib import cm
from matplotlib.colors import Normalize from mpl_toolkits.mplot3d import Axes3D
def sinusoid2d(A, theta, size, T0, Te):
“”” double**2*int*double**2 -> Array[double] “””
ct = np.cos(theta/180*np.pi)
st = np.sin(theta/180*np.pi)
x, y = np.meshgrid(np.arange(0, size, Te), np.arange(0, size, Te)) return A*np.cos(2*np.pi*(y*ct – x*st)/T0)
def shannonInterpolation(I, Te, size):
1

[ ]:
[ ]:
[ ]:
What is the maximal frequency of previous signal s in direction t (denoted fmax) and di- max max max45 t
“”” Array*int*double -> Array[double] “””
n, m = I.shape
x, y = np.meshgrid(np.arange(0, size), np.arange(0, n))
Y = np.sinc(x/Te-y)
x, y = np.meshgrid(np.arange(0, size), np.arange(0, m))
X = np.sinc(x/Te-y)
return np.matmul(X.T, np.matmul(I, Y))
def imshow(I,title=None):
“”” display an image “”” plt.figure(figsize=(500//80,500//80)) plt.gray()
plt.imshow(I)
if title: plt.title(title)
plt.show()
def surf(Z,title=None):
“”” 3D plot of an image “””
X,Y = np.meshgrid(range(Z.shape[1]), range(Z.shape[0]))
fig = plt.figure(figsize=(600/80,600/80))
if title: plt.title(title)
ax = fig.gca(projection=’3d’)
ax.plot_surface(X, Y, Z, cmap=cm.coolwarm, linewidth=0, antialiased=False) plt.show()
1.
2. 2.
2.
We provide the function sinusoid2d(A, theta, L, T0, Te) that allows to sample signal sθ
with a sampling period of Te (the grid is regular with the sample sampling value for direc-
tions u and t). Paremeters A, theta, L and T0 respectively control the amplitude, orientation
and period (T0 = f1 ) of signal sθ. Generate a pseudo continuous signal s45 with A=1, theta 0
= 45,L = 512,T0 = 64andT_e=1.
rection u (denoted fu )? Let fm = max( ft , fu ). Explain why fm is the limit frequency (in sens of Shannon) for the sampling of s45.
(a) Sample s45 with fe = 16 fm and display the sampled signal.
(b) Compute the Fourier transform of the sampled signal and display frequencies. One can use surf() function for an interactive 3D plot.
2

[ ]: [ ]:
[ ]:
[ ]: [ ]:
[ ]:
[ ]:
3.
3.
3.
4.
5.
6.
(b) Write a function error() implementing the relative average error εr = 1LL
2.
(c) Comment the spectrum:
• verify the presence of the two Dirac peaks
• for various values of Te, observe changes in the spectrum. Compare with the spectrum
of the continuous signal (s45). What is the origin of theses differences? • (Bonus question):
– Why, aside the two Dirac peaks, there are somes structures? Explain the origin of these lobes.
– Increase T0 in order to obtain a unique peak. Explain the origin of this fusion. Verify the limit value of T0 for which the two peaks interfer.
(a) Sample s45 with f0 = 4 fm and display the sampled signal.
2AL2 ∑∑|xr(k,l)−xd(k,l)|. k=0 l=0
(c) Reconstruct the sampled signal. Display original and reconstruted signal. Print the relative average error between the original and reconstructed images. What is the origin of this error?
Same question than 3. with fe = 32 fm. Comment the effects of aliasing.
Consider the continuous signal with an oriention of θ = 10. What is the value of fm? With a sampling of 23 fe what is the additional drawback appearing after the reconstruction? Ex- plain.
Bonus question: write a function shannonInterpolationLoop() implementing equation (1) using two loops, in a C way. Compare and discuss the run time of this function and shannonInterpolation() on a small signal (L = 64). Run time way be measured using tic() and tac() functions.
3
from time import process_time
mytime = 0
def tic():
“”” NoneType -> NoneType “””
global mytime

mytime = process_time()
def tac():
“”” NoneType -> int “””
global mytime
print (process_time()-mytime)
mytime = process_time()
### your code starts below
Exercice 2: aliasing on natural images
In this exercice, we study aliasing on image img/barbara.png. Aliasing occurring with subsample of image, we propose to write a code that implements a subsample (using function subSample2() of factor 2 on the image.
[ ]:
from PIL import Image
barbara = np.array(Image.open(‘img/barbara.png’))
def subSample2(I):
“”” Array -> Array “”” return I[::2,::2]
1. Explain what is a subsample of factor 2 and the impact when applied on an image.
2. Write a code that
• iterates the subsampling process
• at each iteration, computes the Fourier transform of the subsampled image
Display subsampled images and their spectrum. Describe and interpret the effects of alias- ing. Why aliasing is a drawback ?
3. Bonus question: same question with the image img/mandrill.png. [ ]:
[ ]:
mandrill = np.array(Image.open(‘img/mandrill.png’)) ### your code and comments start below
4

Practical work 4: Frequency filtering, color
[1]:
[ ]:
[ ]:
2. Write a function idealLowPassFilter(n,m,fc) returning an ideal low pass filter with fre- quency cutoff fc and size n × m. Recall: this function set to 1 pixels at Euclidian distance fc from the center (null frequency).
[ ]:
[ ]:
4. Experiment this function on img/mandrill.png and img/lena.jpg with various values of cut off fc.
• give two effects that appears when fc decreases,
• propose two applications of this filtering.
import numpy as np
from numpy.fft import fft2,ifft2,fftshift
import matplotlib.pyplot as plt
from PIL import Image
def imshow(I,title=None,size=500):
“”” display an image with a specific size “”” plt.figure(figsize=(size//80,size//80)) plt.gray()
plt.imshow(I)
if title: plt.title(title)
plt.show()
Exercice 1 – Frequency filtering
1. Compute than display the centered module of Fourier transform of image img/mandrill.png (use functions seen in previous lessons).
3. Write a function lowPass(I,fc) performing a low pass filtering of an image I. The function should
• compute the centered Fourier transform of I
• multiply point-by-point the spectrum with the ideal low filter produced by
idealLowPassFilter()
• uncenter the filtered spectrum and apply the inverse Fourier transform (function ifft2()
from module numpy.fft
• return the real part of filtered image
1

Exercice 2 – Linear filtering (convolution)
1. Given a kernel convolution of size d × d, d being odd. How many lines and columns should be added to each side of the image to apply this filter? The image is supposed surrounded by zero values.
d/2 lignes d/2 colonnes en plus de chaques côté
2. Write a function imagePad(I,h) getting an image and a kernel, returning a new image padded with zeros according to question 1. It is not allowed to use a module implementing the padding.
[ ]:
[ ]:
4. Try this function on mean filter of size 3 × 3, 5 × 5 and 7 × 7. Discuss the results. [ ]:
5. Display the transfert function of these mean filters. For a better visualization, use the zero- padding technique to obtain a filter with a large size (for instance 256 × 256). Use imshow() and toVisualizeLogFT().
[ ]:
[2]:
3. Write a function conv2(I,h) getting an image and a kernel and returning the convolution of I by h. The function should return an image having the same shape than I. It is not allowed to use a module implementing the convolution.
6. Interpretation: what the analytic expression of the transfert function of a mean filter. Is it an ideal low pass filter?
C’est un sinus cardinal 2d (produit cartésien de deux sinus cardinaux). Non le filtre n’est pas idéal, car à support fréquentiel non borné.
7. Bonus question: perform the same study for the Gaussian kernel. Determine sigma in order to have filter of size 3 × 3, 5 × 5, and 7 × 7.
def gaussianKernel(sigma): “”” double -> Array
return a gaussian kernel of standard deviation sigma
“””
n2 = np.int(3*sigma)
x,y = np.meshgrid(np.arange(-n2,n2+1),np.arange(-n2,n2+1))
kern = np.exp(-(x**2+y**2)/(2*sigma*sigma))
return kern/kern.sum()
### your answer start below
2

Exercice 3: anti aliasing filtering
1. Give a code that subsamples of factor 2 (use function subSample2() given in TME3) the image img/barbara.png.
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
[ ]:
Exercice 4: color image
1. Read images img/clown.bmp and img/clown_lumi.bmp as two arrays named I1 and I2. Dis- play these images examine their shape. What difference there are between them?
2. Give a code that subsamples of factor 2 (use function subSample2()) the image img/barbara.png after applying an low pass filter (use antiAliasingFilter()). As com- ment, recall the principle of filtering in the frequency domain.
def antiAliasingFilter(n,m): “”” int*int -> Array “”” n2, m2 = n//2, m//2
rn, rm = n//4, m//4
A = np.zeros((n, m))
A[rn:rn+n2, rm:rm+m2] = 1
return A
### your answer start below
3. Describe and analyze the filtering of Barbara with and without the anti aliasing filter. What information is lost for the two filtered images ?
2. The first image is an array of dimension 3. Explain the signification of each dimension. From this image create 3 images IR, IG, IB of dimension 2. Display these three images and explain what you see.
3. Create a new image I3 of dimensions 3, the first dimension contains the value of IR, the second the value of IB and the third the values of IG. Try another combinations. Remark: color images are better handled by imshow() if pixel values range in [0, 1].
4. Write a code that allows the see the first channel with red color scales, the second channel in green color scales, and the blue channel in blue color scales.
3

[10]:
Practical work 5: edges detection
The goal of this practial work is to experiment various edges detectors. Attention is given to the following points: 1. comparison between the first and second order detectors 2. study of the impact of smoothing 3. removing non maxima answers of the detectors 4. evaluation in term of robustness and localization
# Useful modules
import numpy as np
from PIL import Image
from scipy.signal import convolve2d
import matplotlib.pyplot as plt
# Useful functions for this work
def orientation(Ix, Iy, Ig):
“”” Array[n,m]**3 -> Array[n,m]
Returns an image of orientation.
“””
n, m = Ix.shape
x = np.arange(4)*np.pi/4
ori = np.stack((np.cos(x), np.sin(x)), axis=1)
O = np.zeros(Ix.shape)
for i in range(n):
for j in range(m):
if Ig[i, j] > 0:
return O
v = np.array([Ix[i, j], -Iy[i, j]])/Ig[i, j]
if Iy[i, j] > 0: v = -v
prod = np.matmul(ori, v)
maxi = prod.max()
imax = np.nonzero(prod == maxi)
O[i, j] = imax[0][0]+1
def gaussianKernel(sigma): “”” double -> Array
return a gaussian kernel of standard deviation sigma
“””
n2 = np.int(3*sigma)
x,y = np.meshgrid(np.arange(-n2,n2+1),np.arange(-n2,n2+1))
kern = np.exp(-(x**2+y**2)/(2*sigma*sigma))
return kern/kern.sum()
def imshow(I, title=None, size=500, axis=False):
“”” display an image, with title, size, and axis “”” plt.figure(figsize=(size//80, size//80))
plt.gray()
plt.imshow(I)
1

if not axis: plt.axis(‘off’)
if title: plt.title(title)
plt.show()
def imshow_hot(I, title=None, size=500, axis=False):
“”” display an image, with title, size, and axis “”” plt.figure(figsize=(size//80, size//80))
plt.hot()
plt.imshow(I)
if not axis: plt.axis(‘off’)
if title: plt.title(title)
plt.show()
def niceDisplay14(affichages,titres=None): “”” list[Array]*list[str] -> NoneType
display from 1 up to 4 images or vectors with optionnal titles
2D arrays are displayed as image with imshow()
1D arrays are displayed as curve with plot()
“””
if not type(affichages) == type([]):
affichages = [affichages]
if titres is None:
titres = [”,]*len(affichages)
if not type(titres) == type([]):
titres = [titres]
nb_affichages = len(affichages)
if nb_affichages >4 or nb_affichages < 1 : raise ValueError('niceDisplay_14 : affichage should be a list of length␣ 􏰀→1 up to 4') if nb_affichages != len(titres): raise ValueError('niceDisplay_14 : titres must have same length than␣ 􏰀→affichage') courbes = False for i in range(0,nb_affichages): s = plt.subplot(101+10*nb_affichages+i) s.set_title(titres[i]) if len(affichages[i].shape)==2 and affichages[i].shape[0] > 1 and␣
􏰀→affichages[i].shape[1] > 1:
# on affiche une image
s.imshow(affichages[i], cmap=”gray”,interpolation=’nearest’,␣ 􏰀→aspect=’equal’)
2

[5]:
[ ]:
else :
# il s’agit d’une seule ligne, à afficher comme une courbe plt.plot(affichages[i])
courbes=True
agrandissement_h = nb_affichages
agrandissement_v = nb_affichages*2 if courbes else nb_affichages
params = plt.gcf()
plSize = params.get_size_inches()
params.set_size_inc