This numerical tour explores some basic image processing tasks.
Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$
#from nt_toolbox.general import *
#from nt_toolbox.signal import *
using NtToolBox
using PyPlot
#using Autoreload # There are some problems with Autoreload package, once solved I will make use of it.
Several functions are implemented to load and display images.
First we load an image.
path to the images
name = "NtToolBox/src/data/lena.png"
n = 256
M = load_image(name, n)
#Mluma = NtToolBox.Mdot(M, [0.299, 0.587, 0.114]) #number of dimensions 2 with the quantity luma.
#Mluma
256×256 Array{Float32,2}: 0.672897 0.672885 0.672751 0.677479 … 0.705561 0.705635 0.528273 0.672897 0.672885 0.672751 0.677479 0.705561 0.705635 0.528273 0.67297 0.672775 0.672727 0.67737 0.704794 0.704138 0.52642 0.668169 0.677114 0.668151 0.663225 0.456002 0.30229 0.195605 0.64486 0.658806 0.663294 0.654061 0.14079 0.126452 0.144403 0.649533 0.649606 0.649768 0.644331 … 0.162321 0.126224 0.153931 0.658769 0.654188 0.654389 0.663061 0.112497 0.135866 0.149514 0.658879 0.663479 0.64497 0.658639 0.134921 0.162588 0.153805 0.649679 0.644896 0.658947 0.65312 0.140041 0.102623 0.168176 0.653712 0.662658 0.658682 0.640679 0.125867 0.143843 0.149978 0.64935 0.644953 0.645078 0.649555 … 0.148765 0.166707 0.162933 0.662948 0.645294 0.635284 0.649588 0.144277 0.180019 0.174482 0.649313 0.653241 0.631026 0.663174 0.162496 0.177566 0.173992 ⋮ ⋱ ⋮ 0.141872 0.241009 0.338395 0.413476 0.274502 0.22575 0.186349 0.121111 0.179189 0.272861 0.33421 … 0.250959 0.152191 0.137347 0.116822 0.154965 0.210891 0.339427 0.175771 0.182736 0.122237 0.136137 0.185606 0.169989 0.239032 0.173045 0.15478 0.180759 0.194887 0.200725 0.200066 0.210942 0.150537 0.195471 0.153815 0.167747 0.197449 0.246806 0.20158 0.141519 0.140765 0.176057 0.149835 0.164028 0.173366 0.150169 … 0.181972 0.181976 0.222797 0.144777 0.173083 0.182443 0.140627 0.213943 0.246564 0.363492 0.135578 0.177404 0.150008 0.163317 0.290465 0.316985 0.373288 0.172575 0.140629 0.16827 0.149863 0.345104 0.36894 0.355407 0.112233 0.149478 0.17255 0.149479 0.38291 0.420395 0.392487 0.116822 0.177807 0.172605 0.149588 … 0.420238 0.42045 0.439224
We can display it. It is possible to zoom on it, extract pixels, etc.
imageplot(M[Int(n/2 - 25) : Int(n/2 + 25), Int(n/2 - 25) : Int(n/2 + 25)], "Zoom", [1, 2, 2])
PyObject <matplotlib.text.Text object at 0x0000000001E4CB38>
An image is a 2D array, that can be modified as a matrix.
imageplot(-M, "-M", [1,2,1])
imageplot(M[end:-1:1,1:size(M, 2)], "Flipped", [1,2,2])
PyObject <matplotlib.text.Text object at 0x000000001E8D6438>
Blurring is achieved by computing a convolution with a kernel.
Compute the low pass Gaussian kernel. Warning, the indexes need to be modulo $n$ in order to use FFTs.
#sigma = 5
#t = concatenate( (arange(0,n/2+1), arange(-n/2,-1)) )
#[Y,X] = np.meshgrid(t,t)
#h = exp( -(X**2+Y**2)/(2.0*float(sigma)**2) )
#h = h/sum(h)
#imageplot(fftshift(h))
sigma = 5
X = [0:n/2; -n/2:-2]'
Y = [0:n/2; -n/2:-2]
h = exp(-(X.^2 .+ Y.^2)/(2*(sigma)^2))
h = h/sum(h)
imageplot(fftshift(h))
h
256×256 Array{Float64,2}: 0.00749229 0.00734394 0.00691626 … 0.00625809 0.00691626 0.00734394 0.00719852 0.00677931 0.00613417 0.00677931 0.00691626 0.00677931 0.00638451 0.00577694 0.00638451 0.00625809 0.00613417 0.00577694 0.0052272 0.00577694 0.00544052 0.00533279 0.00502223 0.00454431 0.00502223 0.00454431 0.00445432 0.00419492 … 0.00379572 0.00419492 0.00364689 0.00357468 0.0033665 0.00304614 0.0033665 0.00281194 0.00275626 0.00259575 0.00234873 0.00259575 0.00208314 0.00204189 0.00192298 0.00173998 0.00192298 0.00148272 0.00145336 0.00136872 0.00123847 0.00136872 0.00101397 0.000993894 0.000936014 … 0.00084694 0.000936014 0.000666227 0.000653035 0.000615005 0.000556479 0.000615005 0.000420578 0.00041225 0.000388243 0.000351296 0.000388243 ⋮ ⋱ ⋮ 0.000255094 0.000250042 0.000235481 0.000213072 0.000235481 0.000420578 0.00041225 0.000388243 … 0.000351296 0.000388243 0.000666227 0.000653035 0.000615005 0.000556479 0.000615005 0.00101397 0.000993894 0.000936014 0.00084694 0.000936014 0.00148272 0.00145336 0.00136872 0.00123847 0.00136872 0.00208314 0.00204189 0.00192298 0.00173998 0.00192298 0.00281194 0.00275626 0.00259575 … 0.00234873 0.00259575 0.00364689 0.00357468 0.0033665 0.00304614 0.0033665 0.00454431 0.00445432 0.00419492 0.00379572 0.00419492 0.00544052 0.00533279 0.00502223 0.00454431 0.00502223 0.00625809 0.00613417 0.00577694 0.0052272 0.00577694 0.00691626 0.00677931 0.00638451 … 0.00577694 0.00638451
Compute the periodic convolution ussing FFTs
#Mh = real( plan_ifft(plan_fft(M) * plan_fft(h)) )
#Mh = conv(M,h)
#Mh1 = conv2(Array{Float64, 2}(M[:, :, 1]), h)
#Mh2 = conv2(Array{Float64, 2}(M[:, :, 2]), h)
#Mh3 = conv2(Array{Float64, 2}(M[:, :, 3]), h)
#Mh1 = Mh1[1:256, 1:256]
#Mh2 = Mh2[1:256, 1:256]
#Mh3 = Mh3[1:256, 1:256]
#Mh = zeros(n, n, 3)
#Mh[:, :, 1] = Mh1
#Mh[:, :, 2] = Mh2
#Mh[:, :, 3] = Mh3
Mh = conv2(Array{Float64, 2}(M), h)
#Mh = Mh[1:256, 1:256]
Mh = Mh[1:255, 1:255] + Mh[257:511, 1:255] + Mh[1:255, 257:511] + Mh[257:511, 257:511]
255×255 Array{Float64,2}: 0.399904 0.414308 0.430863 0.449494 … 0.365755 0.370753 0.378107 0.419226 0.434839 0.452145 0.470926 0.380194 0.386088 0.394598 0.434686 0.451881 0.470352 0.489711 0.39024 0.396975 0.406717 0.445423 0.464554 0.484593 0.504965 0.395006 0.402559 0.413623 0.451573 0.47292 0.494845 0.516585 0.394526 0.402973 0.415504 0.453679 0.477411 0.501421 0.524766 … 0.389342 0.398844 0.413014 0.452427 0.47863 0.504826 0.529918 0.380193 0.39096 0.406945 0.448729 0.477382 0.505757 0.532627 0.368138 0.3804 0.398355 0.443463 0.474477 0.504951 0.533556 0.354215 0.368191 0.388236 0.437565 0.470776 0.50319 0.533407 0.339613 0.355477 0.377674 0.431751 0.466934 0.501074 0.532728 … 0.325326 0.343182 0.367514 0.426451 0.463369 0.499009 0.53191 0.311965 0.331849 0.358248 0.421951 0.460362 0.497278 0.531235 0.30002 0.321887 0.35022 ⋮ ⋱ 0.295707 0.311493 0.332376 0.357719 0.284313 0.280207 0.280355 0.292591 0.30854 0.329552 0.35502 0.279993 0.276231 0.276727 0.288856 0.304762 0.325715 0.351162 … 0.275598 0.272112 0.272838 0.285443 0.301132 0.321846 0.347111 0.27181 0.268606 0.269516 0.283304 0.298629 0.318932 0.343837 0.26938 0.266526 0.267625 0.283363 0.298234 0.317995 0.342378 0.268987 0.266602 0.26796 0.286396 0.30077 0.319893 0.343608 0.271267 0.269495 0.271208 0.293053 0.306928 0.325351 0.348267 … 0.276766 0.275757 0.277944 0.303685 0.31712 0.334832 0.356852 0.28575 0.285643 0.288431 0.318179 0.331298 0.348356 0.369445 0.298072 0.298968 0.302481 0.336006 0.348998 0.365533 0.385728 0.313213 0.315161 0.319502 0.356286 0.369411 0.385638 0.405066 0.330278 0.333281 0.33854
Display
imageplot(M, "Image", [1, 2, 1])
imageplot(Mh, "Blurred", [1, 2, 2])
PyObject <matplotlib.text.Text object at 0x000000001E7DCDA0>
Several differential and convolution operators are implemented.
(G_x, G_y) = Images.imgradients(M)
imageplot(G_x, "d/ dx", [1, 2, 1])
imageplot(G_y, "d/ dy", [1, 2, 2])
WARNING: the order of outputs has switched (`grad1, grad2 = imgradients(img)` rather than `gradx, grady = imgradients`). Silence this warning by providing a kernelfun, e.g., imgradients(img, KernelFactors.ando3).
in depwarn(::String, ::Symbol) at .\deprecated.jl:64
in imgradients(::Array{Float32,2}) at C:\Users\Ayman\.julia\v0.5\ImageFiltering\src\specialty.jl:50
in include_string(::String, ::String) at .\loading.jl:441
in execute_request(::ZMQ.Socket, ::IJulia.Msg) at C:\Users\Ayman\.julia\v0.5\IJulia\src\execute_request.jl:157
in eventloop(::ZMQ.Socket) at C:\Users\Ayman\.julia\v0.5\IJulia\src\eventloop.jl:8
in
(::IJulia.##13#19)() at .\task.jl:360 while loading In[8], in expression starting on line 1
PyObject <matplotlib.text.Text object at 0x00000000228428D0>
The 2D Fourier transform can be used to perform low pass approximation and interpolation (by zero padding).
Compute and display the Fourier transform (display over a log scale). The function fftshift is useful to put the 0 low frequency in the middle. After fftshift, the zero frequency is located at position $(n/2+1,n/2+1)$.
Mf = plan_fft(M)
Mf*M #To verify.
Lf = fftshift(log(abs(Mf*M) + 1e-1))
imageplot(M, "Image", [1, 2, 1])
imageplot(Lf, "Fourier transform", [1, 2, 2])
PyObject <matplotlib.text.Text object at 0x0000000022C59128>
Exercise 1: To avoid boundary artifacts and estimate really the frequency content of the image (and not of the artifacts!), one needs to multiply M by a smooth windowing function h and compute fft2(M*h). Use a sine windowing function. Can you interpret the resulting filter ?
#run -i nt_solutions/introduction_3_image/exo1
include("introduction_3_image\ exo1.jl")
PyObject <matplotlib.text.Text object at 0x0000000029589668>
Exercise 2: Perform low pass filtering by removing the high frequencies of the spectrum. What do you oberve ?
#run -i nt_solutions/introduction_3_image/exo2
include("introduction_3_image\ exo2.jl")
PyObject <matplotlib.text.Text object at 0x0000000029940DD8>