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}$
This numerical tour overviews the use of Fourier and wavelets for image approximation.
#from __future__ import division
#import numpy as np
#import scipy as scp
#import pylab as pyl
#import matplotlib.pyplot as plt
#from nt_toolbox.general import *
#from nt_toolbox.signal import *
#import warnings
#warnings.filterwarnings('ignore')
#%matplotlib inline
#%load_ext autoreload
#%autoreload 2
using PyPlot
using NtToolBox
using Autoreload
arequire("NtToolBox")
WARNING: symbol is deprecated, use Symbol instead. in depwarn(::String, ::Symbol) at .\deprecated.jl:64 in symbol(::String, ::Vararg{String,N}) at .\deprecated.jl:30 in reload_module(::Symbol, ::Expr) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\smart_types.jl:65 in (::Autoreload.##14#15{Array{Any,1}})() at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:117 in cd(::Autoreload.##14#15{Array{Any,1}}, ::String) at .\file.jl:48 in #smart_reload#13(::Array{Any,1}, ::Function, ::String) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:109 in (::Autoreload.#kw##smart_reload)(::Array{Any,1}, ::Autoreload.#smart_reload, ::String) at .\<missing>:0 in #try_reload#16(::Array{Any,1}, ::Function, ::String) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:126 in (::Autoreload.#kw##try_reload)(::Array{Any,1}, ::Autoreload.#try_reload, ::String) at .\<missing>:0 in #areload#17(::Array{Any,1}, ::Function, ::Symbol) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:185 in areload_hook() at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:197 in execute_request(::ZMQ.Socket, ::IJulia.Msg) at C:\Users\Ayman\.julia\v0.5\IJulia\src\execute_request.jl:153 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 C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl, in expression starting on line 31 WARNING: Main.(:NtToolBox) is deprecated; use Main.:NtToolBox or getfield(Main, :NtToolBox) instead. in depwarn(::String, ::Symbol) at .\deprecated.jl:64 in broadcast(::Module, ::Symbol) at .\deprecated.jl:377 in reload_module(::Symbol, ::Expr) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\smart_types.jl:72 in (::Autoreload.##14#15{Array{Any,1}})() at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:117 in cd(::Autoreload.##14#15{Array{Any,1}}, ::String) at .\file.jl:48 in #smart_reload#13(::Array{Any,1}, ::Function, ::String) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:109 in (::Autoreload.#kw##smart_reload)(::Array{Any,1}, ::Autoreload.#smart_reload, ::String) at .\<missing>:0 in #try_reload#16(::Array{Any,1}, ::Function, ::String) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:126 in (::Autoreload.#kw##try_reload)(::Array{Any,1}, ::Autoreload.#try_reload, ::String) at .\<missing>:0 in #areload#17(::Array{Any,1}, ::Function, ::Symbol) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:185 in areload_hook() at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:197 in execute_request(::ZMQ.Socket, ::IJulia.Msg) at C:\Users\Ayman\.julia\v0.5\IJulia\src\execute_request.jl:153 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 C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl, in expression starting on line 31 WARNING: Method definition Mdot
perform_wavelet_transf (generic function with 4 methods)
(Any, Array{Float64, 1}) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:15 overwritten at none:16. WARNING: Method definition rescale(Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:29 overwritten at none:30. WARNING: Method definition rescale(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:29 overwritten at none:30. WARNING: Method definition rescale(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:29 overwritten at none:30. WARNING: Method definition clamP(Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:43 overwritten at none:44. WARNING: Method definition clamP(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:43 overwritten at none:44. WARNING: Method definition clamP(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:43 overwritten at none:44. WARNING: Method definition gaussian_blur(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:78 overwritten at none:79. WARNING: Method definition plot_levelset(Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:99 overwritten at none:100. WARNING: Method definition plot_levelset(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:99 overwritten at none:100. WARNING: Method definition plot_levelset(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:99 overwritten at none:100. WARNING: Method definition Grad(Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:113 overwritten at none:114. WARNING: Method definition Grad(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:113 overwritten at none:114. WARNING: Method definition Grad(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\NtToolBox.jl:113 overwritten at none:114. WARNING: Method definition ndgrid(AbstractArray{T<:Any, 1}) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:3 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:3. WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:6 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:6. WARNING: Method definition ndgrid_fill(Any, Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:13 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:13. WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}...) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:19 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:19. WARNING: Method definition meshgrid(AbstractArray{T<:Any, 1}) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:33 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:33. WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:36 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:36. WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:44 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\ndgrid.jl:44. WARNING: Method definition load_image(Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:12 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:12. WARNING: Method definition load_image(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:12 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:12. WARNING: Method definition load_image(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:12 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:12. WARNING: Method definition load_image(Any, Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:12 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:12. WARNING: Method definition load_image(Any, Any, Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:12 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:12. WARNING: Method definition imageplot(Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:51 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:51. WARNING: Method definition imageplot(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:51 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:51. WARNING: Method definition imageplot(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:51 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:51. WARNING: Method definition snr(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:63 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:63. WARNING: Method definition resize_image(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:81 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:81. WARNING: Method definition trim_dim(Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:100 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:100. WARNING: Method definition subsampling(Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:120 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:120. WARNING: Method definition subsampling(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:120 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:120. WARNING: Method definition subsampling(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:120 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:120. WARNING: Method definition upsampling(Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:146 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:146. WARNING: Method definition upsampling(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:146 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:146. WARNING: Method definition upsampling(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:146 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:146. WARNING: Method definition cconvol(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:178 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:178. WARNING: Method definition cconvol(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:178 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:178. WARNING: Method definition plot_wavelet(Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:255 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:255. WARNING: Method definition plot_wavelet(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:255 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:255. WARNING: Method definition subselectdim(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:306 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:306. WARNING: Method definition subassign(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:320 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:320. WARNING: Method definition subselect(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:334 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:334. WARNING: Method definition perform_wavortho_transf(Any, Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:359 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:359. WARNING: Method definition div(Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:394 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:394. WARNING: Method definition grad(Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:406 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:406. WARNING: Method definition bilinear_interpolate(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:418 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\signal.jl:418. WARNING: Method definition perform_wavelet_transf(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:3 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:3. WARNING: Method definition perform_wavelet_transf(Any, Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:3 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:3. WARNING: Method definition perform_wavelet_transf(Any, Any, Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:3 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:3. WARNING: Method definition perform_wavelet_transf(Any, Any, Any, Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:3 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:3. WARNING: Method definition lifting_step(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:158 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:158. WARNING: Method definition lifting_step_ti(Any, Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:200 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:200. WARNING: Method definition compute_wavelet_filter(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\compute_wavelet_filter.jl:2 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\compute_wavelet_filter.jl:2. WARNING: Method definition Div(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\Div.jl:2 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\Div.jl:2. WARNING: Method definition Div(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\Div.jl:2 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\Div.jl:2. WARNING: Method definition Div(Any, Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\Div.jl:2 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\Div.jl:2. WARNING: Method definition perform_blurring(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:3 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:3. WARNING: Method definition perform_blurring(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:3 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:3. WARNING: Method definition compute_gaussian_filter(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:43 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:43. WARNING: Method definition build_gaussian_filter_2d(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:85 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:85. WARNING: Method definition build_gaussian_filter_2d(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:85 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:85. WARNING: Method definition build_gaussian_filter_1d(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:132 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:132. WARNING: Method definition build_gaussian_filter_1d(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:132 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:132. WARNING: Method definition perform_convolution(Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:161 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:161. WARNING: Method definition perform_convolution(Any, Any, Any) in module NtToolBox at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:161 overwritten at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_blurring.jl:161. WARNING: using NtToolBox.plot_wavelet in module Main conflicts with an existing identifier. WARNING: using NtToolBox.imageplot in module Main conflicts with an existing identifier. WARNING: using NtToolBox.snr in module Main conflicts with an existing identifier. WARNING: using NtToolBox.load_image in module Main conflicts with an existing identifier. WARNING: using NtToolBox.rescale in module Main conflicts with an existing identifier. WARNING: using NtToolBox.perform_wavelet_transf in module Main conflicts with an existing identifier. WARNING: using NtToolBox.clamP in module Main conflicts with an existing identifier. WARNING: replacing module NtToolBox
Note: to measure the error of an image $f$ with its approximation $f_M$, we use the SNR measure, defined as
$$ \text{SNR}(f,f_M) = -20\log_{10} \pa{ \frac{ \norm{f-f_M} }{ \norm{f} } }, $$which is a quantity expressed in decibels (dB). The higer the SNR, the better the quality.
First we load an image $ f \in \RR^N $ of $ N = N_0 \times N_0 $ pixels.
n0 = 512
f = rescale(load_image("NtToolBox/src/data/lena.png", n0))
512×512 Array{Float32,2}: 0.672897 0.672897 0.672897 0.669782 … 0.705608 0.64486 0.528037 0.672897 0.672897 0.672897 0.669782 0.705608 0.64486 0.528037 0.672897 0.672897 0.672897 0.669782 0.705608 0.64486 0.528037 0.672897 0.672897 0.672897 0.669782 0.705608 0.64486 0.528037 0.672897 0.672897 0.672897 0.669782 0.705608 0.64486 0.528037 0.682243 0.682243 0.658879 0.64486 … 0.514019 0.415888 0.291277 0.668224 0.668224 0.67757 0.658879 0.303738 0.224299 0.196262 0.663551 0.663551 0.64486 0.654206 0.149533 0.17757 0.135514 0.64486 0.64486 0.658879 0.658879 0.126168 0.135514 0.14486 0.64486 0.64486 0.654206 0.658879 0.135514 0.126168 0.116822 0.649533 0.649533 0.649533 0.668224 … 0.126168 0.154206 0.154206 0.649533 0.649533 0.649533 0.663551 0.126168 0.121495 0.140187 0.658879 0.658879 0.654206 0.649533 0.135514 0.140187 0.149533 ⋮ ⋱ ⋮ 0.163551 0.163551 0.182243 0.200935 … 0.17757 0.168224 0.154206 0.149533 0.149533 0.163551 0.182243 0.182243 0.228972 0.224299 0.140187 0.140187 0.182243 0.219626 0.200935 0.257009 0.317757 0.14486 0.14486 0.172897 0.17757 0.247664 0.280374 0.364486 0.140187 0.140187 0.168224 0.14486 0.261682 0.299065 0.336449 0.135514 0.135514 0.17757 0.168224 … 0.317757 0.359813 0.373832 0.140187 0.140187 0.186916 0.163551 0.35514 0.345794 0.378505 0.172897 0.172897 0.140187 0.135514 0.369159 0.373832 0.35514 0.126168 0.126168 0.149533 0.149533 0.411215 0.392523 0.383178 0.11215 0.11215 0.149533 0.135514 0.420561 0.401869 0.392523 0.116822 0.116822 0.17757 0.154206 … 0.420561 0.425234 0.439252 0.116822 0.116822 0.17757 0.238318 0.420561 0.425234 0.439252
Display the original image.
figure(figsize = (5,5))
imageplot(f, "Image f")
PyObject <matplotlib.text.Text object at 0x000000001E519B00>
Display a zoom in the middle.
figure(figsize = (5,5))
imageplot(f[Int(n0/2 - 32) : Int(n0/2 + 32), Int(n0/2 - 32) : Int(n0/2 + 32)], "Zoom")
PyObject <matplotlib.text.Text object at 0x000000001E574C88>
An image is a 2D array, it can be modified as a matrix.
figure(figsize = (8,8))
imageplot(-f, "-f", [1, 2, 1])
imageplot(f[end:-1:1, :], "Flipped", [1, 2, 2])
PyObject <matplotlib.text.Text object at 0x000000001EF9B400>
Blurring is achieved by computing a convolution $f \star h$ with a kernel $h$.
Compute the low pass kernel.
k = 9; #size of the kernel
h = ones(k, k)
h = h/sum(h) #normalize
9×9 Array{Float64,2}: 0.0123457 0.0123457 0.0123457 … 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 … 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457 0.0123457
Compute the convolution $f \star h$.
fh = conv2(Array{Float64, 2}(f), h)
520×520 Array{Float64,2}: 0.00830737 0.0166147 0.0249221 … 0.0231914 0.0144802 0.00651898 0.0166147 0.0332295 0.0498442 0.0463828 0.0289604 0.013038 0.0249221 0.0498442 0.0747664 0.0695743 0.0434406 0.0195569 0.0332295 0.066459 0.0996885 0.0927657 0.0579209 0.0260759 0.0415369 0.0830737 0.124611 0.115957 0.0724011 0.0325949 0.0499596 0.0999192 0.14959 … 0.131033 0.0811315 0.0361909 0.0582093 0.116419 0.174455 0.139975 0.0863236 0.0386139 0.0664013 0.132803 0.1988 0.145687 0.0901888 0.0402869 0.0743625 0.148725 0.222857 0.150706 0.0936502 0.0420753 0.0740164 0.148033 0.221934 0.132187 0.0821699 0.0369986 0.0737279 0.147456 0.221068 … 0.114361 0.0714973 0.0323834 0.0734395 0.146879 0.220203 0.0959578 0.0602477 0.0275951 0.0732664 0.146533 0.219626 0.0780162 0.0493443 0.0229222 ⋮ ⋱ 0.0162109 0.0324218 0.0512288 0.100208 0.069055 0.0356525 0.0155763 0.0311526 0.0495558 0.109034 0.0748817 0.0385947 0.0151725 0.030345 0.0489212 … 0.117053 0.0799585 0.0412484 0.014884 0.0297681 0.0482866 0.123341 0.0835352 0.0427484 0.0130956 0.0261913 0.0425753 0.112323 0.075574 0.0382485 0.0113649 0.0227299 0.037037 0.101246 0.0677282 0.0340948 0.00969193 0.0193839 0.0314988 0.0882658 0.0586708 0.0294796 0.00796123 0.0159225 0.0257298 … 0.0749394 0.0497289 0.0248067 0.0058267 0.0116534 0.01973 0.0613823 0.0407292 0.0204223 0.00426907 0.00853813 0.0147687 0.046729 0.0311526 0.0156917 0.0028845 0.00576901 0.0101535 0.0317295 0.0213453 0.0108457 0.00144225 0.0028845 0.00507673 0.0158648 0.0106727 0.00542287
Display.
figure(figsize = (5,5))
imageplot(fh, "Blurred image")
PyObject <matplotlib.text.Text object at 0x000000001F088E48>
The Fourier orthonormal basis is defined as $$ \psi_m(k) = \frac{1}{\sqrt{N}}e^{\frac{2i\pi}{N_0} \dotp{m}{k} } $$ where $0 \leq k_1,k_2 < N_0$ are position indexes, and $0 \leq m_1,m_2 < N_0$ are frequency indexes.
The Fourier transform $\hat f$ is the projection of the image on this Fourier basis
$$ \hat f(m) = \dotp{f}{\psi_m}. $$The Fourier transform is computed in $ O(N \log(N)) $ operation using the FFT algorithm (Fast Fourier Transform). Note the normalization by $\sqrt{N}=N_0$ to make the transform orthonormal.
F = plan_fft(f)
F = (F*f)/n0
512×512 Array{Complex{Float32},2}: 256.476+0.0im -3.95615+20.4749im … -3.95615-20.4749im 0.0844855-11.1962im -15.1921+13.9858im 0.630981+7.94724im -3.55624-1.76501im -2.36528-0.595752im 7.28666-0.598839im 4.73903+1.55932im -0.805839-5.27302im 6.26791-6.04117im -0.624841-5.30147im 1.55698-1.17529im -0.988565+0.720974im -1.61084-4.77705im 0.165628-1.66687im … -2.28491+1.02422im 1.22427-0.125188im -2.10585-2.96785im 1.15793-2.06182im 2.86644-0.98837im -0.248534+1.02317im -0.55987-1.09585im -0.569311-0.075523im 0.416833-0.219624im 1.58716-0.4015im 0.993204-1.84695im -1.48726-0.244343im -1.33509+0.742391im -0.0138412-0.112531im -0.682093-1.73745im … 1.08995-0.3994im 0.83187-1.05805im 0.247246+0.45962im 0.700845-0.604206im -0.225175+0.184748im -0.398911-0.745511im 1.23564-0.586102im ⋮ ⋱ -0.225175-0.184748im 1.23564+0.586102im … -0.398911+0.745511im 0.83187+1.05805im 0.700845+0.604206im 0.247246-0.45962im -0.0138412+0.112531im 1.08995+0.3994im -0.682093+1.73745im 0.993204+1.84695im -1.33509-0.742392im -1.48726+0.244343im -0.569311+0.075523im 1.58716+0.4015im 0.416833+0.219624im 2.86644+0.98837im -0.55987+1.09585im … -0.248534-1.02317im 1.22427+0.125188im 1.15793+2.06182im -2.10585+2.96785im -1.61084+4.77705im -2.28491-1.02422im 0.165628+1.66687im -0.624841+5.30147im -0.988565-0.720973im 1.55698+1.17529im 4.73903-1.55932im 6.26791+6.04117im -0.805839+5.27302im -3.55624+1.76501im 7.28667+0.598839im … -2.36528+0.595752im 0.0844855+11.1962im 0.630981-7.94724im -15.1921-13.9858im
We check this conservation of the energy.
println(@sprintf("Energy of Image: %f", norm(f)))
println(@sprintf("Energy of Fourier: %f", norm(F)))
Energy of Image: 262.554108 Energy of Fourier: 262.554138
Compute the logarithm of the Fourier magnitude $ \log\left(\abs{\hat f(m)} + \epsilon\right) $, for some small $\epsilon$.
L = fftshift(log(abs(F) + 1e-1))
512×512 Array{Float32,2}: -2.22263 -2.13613 -2.2846 -2.2275 … -2.2275 -2.2846 -2.13613 -2.20524 -2.25339 -2.26194 -2.14533 -2.2274 -2.20701 -2.16727 -2.20634 -2.23409 -2.19394 -2.29619 -2.16088 -2.25066 -2.24465 -2.21319 -2.28514 -2.25704 -2.25984 -2.11008 -2.19926 -2.19242 -2.21663 -2.20098 -2.2657 -2.20428 -2.21889 -2.17542 -2.2941 -2.22902 -2.17697 -2.22038 -2.19271 … -2.28654 -2.17516 -2.13996 -2.22738 -2.1867 -2.14879 -2.20321 -2.26124 -2.20302 -2.11106 -2.24976 -2.19271 -2.22961 -2.21817 -2.17068 -2.25378 -2.22941 -2.16005 -2.11057 -2.27793 -2.08118 -2.22947 -2.23754 -2.28634 -2.1473 -2.18849 -2.18699 -2.15036 -2.1783 -2.29059 -2.2841 -2.25382 -2.16019 -2.22127 -2.19619 … -2.27501 -2.20996 -2.2211 -2.24757 -2.20535 -2.25617 -2.21006 -2.19404 -2.09374 -2.17059 -2.17976 -2.17071 -2.13844 -2.23593 -2.26239 -2.25548 -2.21712 ⋮ ⋱ ⋮ -2.17976 -2.21712 -2.25548 -2.26239 … -2.23593 -2.13844 -2.17071 -2.24757 -2.17059 -2.09374 -2.19404 -2.21006 -2.25617 -2.20535 -2.25382 -2.2211 -2.20996 -2.27501 -2.19619 -2.22127 -2.16019 -2.1473 -2.2841 -2.29059 -2.1783 -2.15036 -2.18699 -2.18849 -2.16005 -2.28634 -2.23754 -2.22947 -2.08118 -2.27793 -2.11057 -2.24976 -2.22941 -2.25378 -2.17068 … -2.21817 -2.22961 -2.19271 -2.22738 -2.11106 -2.20302 -2.26124 -2.20321 -2.14879 -2.1867 -2.22902 -2.13996 -2.17516 -2.28654 -2.19271 -2.22038 -2.17697 -2.21663 -2.2941 -2.17542 -2.21889 -2.20428 -2.2657 -2.20098 -2.21319 -2.19242 -2.19926 -2.11008 -2.25984 -2.25704 -2.28514 -2.20634 -2.24465 -2.25066 -2.16088 … -2.29619 -2.19394 -2.23409 -2.20524 -2.16727 -2.20701 -2.2274 -2.14532 -2.26194 -2.25339
Display. Note that we use the function fftshift to put the 0 low frequency in the middle.
figure(figsize = (5,5))
imageplot(L, "Log(Fourier transform)")
PyObject <matplotlib.text.Text object at 0x000000002805DB00>
An approximation is obtained by retaining a certain set of index $I_M$
$$ f_M = \sum_{ m \in I_M } \dotp{f}{\psi_m} \psi_m. $$Linear approximation is obtained by retaining a fixed set $I_M$ of $M = \abs{I_M}$ coefficients. The important point is that $I_M$ does not depend on the image $f$ to be approximated.
For the Fourier transform, a low pass linear approximation is obtained by keeping only the frequencies within a square.
$$ I_M = \enscond{m=(m_1,m_2)}{ -q/2 \leq m_1,m_2 < q/2 } $$where $ q = \sqrt{M} $.
This can be achieved by computing the Fourier transform, setting to zero the $N-M$ coefficients outside the square $I_M$ and then inverting the Fourier transform.
Number $M$ of kept coefficients.
M = Base.div(n0^2, 64)
4096
Exercise 1
Perform the linear Fourier approximation with $M$ coefficients. Store the result in the variable $f_M$.
#run -i nt_solutions/introduction_4_fourier_wavelets/exo1
include("introduction_4_fourier_wavelets\ exo1.jl")
PyObject <matplotlib.text.Text object at 0x0000000028328588>
## Insert your code here.
Compare two 1D profile (lines of the image). This shows the strong ringing artifact of the linea approximation.
figure(figsize = (7, 6))
subplot(2, 1, 1)
plot(f[: , Base.div(n0, 2)])
xlim(0, n0)
title("f")
subplot(2, 1, 2)
plot(fM[: , Base.div(n0, 2)])
xlim(0, n0)
title("f_M")
show()
Non-linear approximation is obtained by keeping the $M$ largest coefficients. This is equivalently computed using a thresholding of the coefficients $$ I_M = \enscond{m}{ \abs{\dotp{f}{\psi_m}}>T }. $$
Set a threshold $T>0$.
T = .2
0.2
Compute the Fourier transform.
F = plan_fft(f)
F = (F*f)/n0
512×512 Array{Complex{Float32},2}: 256.476+0.0im -3.95615+20.4749im … -3.95615-20.4749im 0.0844855-11.1962im -15.1921+13.9858im 0.630981+7.94724im -3.55624-1.76501im -2.36528-0.595752im 7.28666-0.598839im 4.73903+1.55932im -0.805839-5.27302im 6.26791-6.04117im -0.624841-5.30147im 1.55698-1.17529im -0.988565+0.720974im -1.61084-4.77705im 0.165628-1.66687im … -2.28491+1.02422im 1.22427-0.125188im -2.10585-2.96785im 1.15793-2.06182im 2.86644-0.98837im -0.248534+1.02317im -0.55987-1.09585im -0.569311-0.075523im 0.416833-0.219624im 1.58716-0.4015im 0.993204-1.84695im -1.48726-0.244343im -1.33509+0.742391im -0.0138412-0.112531im -0.682093-1.73745im … 1.08995-0.3994im 0.83187-1.05805im 0.247246+0.45962im 0.700845-0.604206im -0.225175+0.184748im -0.398911-0.745511im 1.23564-0.586102im ⋮ ⋱ -0.225175-0.184748im 1.23564+0.586102im … -0.398911+0.745511im 0.83187+1.05805im 0.700845+0.604206im 0.247246-0.45962im -0.0138412+0.112531im 1.08995+0.3994im -0.682093+1.73745im 0.993204+1.84695im -1.33509-0.742392im -1.48726+0.244343im -0.569311+0.075523im 1.58716+0.4015im 0.416833+0.219624im 2.86644+0.98837im -0.55987+1.09585im … -0.248534-1.02317im 1.22427+0.125188im 1.15793+2.06182im -2.10585+2.96785im -1.61084+4.77705im -2.28491-1.02422im 0.165628+1.66687im -0.624841+5.30147im -0.988565-0.720973im 1.55698+1.17529im 4.73903-1.55932im 6.26791+6.04117im -0.805839+5.27302im -3.55624+1.76501im 7.28667+0.598839im … -2.36528+0.595752im 0.0844855+11.1962im 0.630981-7.94724im -15.1921-13.9858im
Do the hard thresholding.
#abs(F) .> T
FT = F .* (abs(F) .> T)
512×512 Array{Complex{Float32},2}: 256.476+0.0im -3.95615+20.4749im … -3.95615-20.4749im 0.0844855-11.1962im -15.1921+13.9858im 0.630981+7.94724im -3.55624-1.76501im -2.36528-0.595752im 7.28666-0.598839im 4.73903+1.55932im -0.805839-5.27302im 6.26791-6.04117im -0.624841-5.30147im 1.55698-1.17529im -0.988565+0.720974im -1.61084-4.77705im 0.165628-1.66687im … -2.28491+1.02422im 1.22427-0.125188im -2.10585-2.96785im 1.15793-2.06182im 2.86644-0.98837im -0.248534+1.02317im -0.55987-1.09585im -0.569311-0.075523im 0.416833-0.219624im 1.58716-0.4015im 0.993204-1.84695im -1.48726-0.244343im -1.33509+0.742391im -0.0-0.0im -0.682093-1.73745im … 1.08995-0.3994im 0.83187-1.05805im 0.247246+0.45962im 0.700845-0.604206im -0.225175+0.184748im -0.398911-0.745511im 1.23564-0.586102im ⋮ ⋱ -0.225175-0.184748im 1.23564+0.586102im … -0.398911+0.745511im 0.83187+1.05805im 0.700845+0.604206im 0.247246-0.45962im -0.0+0.0im 1.08995+0.3994im -0.682093+1.73745im 0.993204+1.84695im -1.33509-0.742392im -1.48726+0.244343im -0.569311+0.075523im 1.58716+0.4015im 0.416833+0.219624im 2.86644+0.98837im -0.55987+1.09585im … -0.248534-1.02317im 1.22427+0.125188im 1.15793+2.06182im -2.10585+2.96785im -1.61084+4.77705im -2.28491-1.02422im 0.165628+1.66687im -0.624841+5.30147im -0.988565-0.720973im 1.55698+1.17529im 4.73903-1.55932im 6.26791+6.04117im -0.805839+5.27302im -3.55624+1.76501im 7.28667+0.598839im … -2.36528+0.595752im 0.0844855+11.1962im 0.630981-7.94724im -15.1921-13.9858im
Display. Note that we use the function fftshift to put the 0 low frequency in the middle.
L = fftshift(log(abs(FT) + 1e-1))
figure(figsize = (5,5))
imageplot(L, "thresholded Log(Fourier transform)")
PyObject <matplotlib.text.Text object at 0x00000000289AB358>
Inverse Fourier transform to obtain $f_M$.
fM = plan_ifft(FT)
fM = real(n0*(fM*FT))
512×512 Array{Float32,2}: 0.487353 0.461108 0.442373 0.436895 … 0.585559 0.548152 0.515611 0.54188 0.529753 0.520657 0.518191 0.625147 0.585797 0.558098 0.577563 0.580702 0.582367 0.584125 0.637712 0.600003 0.579998 0.591426 0.61124 0.624315 0.631086 0.6142 0.58382 0.576778 0.585582 0.623212 0.647824 0.660066 0.554214 0.537927 0.550134 0.565794 0.621657 0.657638 0.675643 … 0.467046 0.470584 0.506988 0.539474 0.61319 0.660248 0.683929 0.368853 0.394732 0.456941 0.513838 0.604336 0.661789 0.689986 0.277286 0.323711 0.409414 0.494351 0.59982 0.665919 0.695777 0.205958 0.267419 0.371379 0.483499 0.601162 0.672627 0.699845 0.160904 0.230096 0.345998 0.480211 0.606355 0.678847 0.699168 … 0.14001 0.210078 0.332255 0.480526 0.61123 0.680914 0.692205 0.13511 0.201267 0.325739 0.479675 0.612156 0.67754 0.681115 0.135659 0.195813 0.320626 ⋮ ⋱ ⋮ 0.206354 0.202004 0.195751 0.189541 … 0.160972 0.183008 0.200854 0.208493 0.198254 0.19192 0.187799 0.190108 0.207397 0.214143 0.209372 0.193383 0.187316 0.186153 0.221985 0.231609 0.226227 0.210642 0.187599 0.179899 0.180694 0.253315 0.255072 0.238493 0.21288 0.180328 0.167683 0.16842 0.282435 0.277903 0.251839 0.216362 0.17228 0.151818 0.150757 … 0.309282 0.300268 0.266377 0.222645 0.167318 0.13829 0.13479 0.334999 0.322428 0.282303 0.235526 0.172405 0.136685 0.131131 0.361742 0.345515 0.301191 0.260183 0.195097 0.156346 0.149512 0.392584 0.372167 0.32636 0.300544 0.239779 0.202002 0.194524 0.430826 0.405866 0.361457 0.356369 0.304823 0.271296 0.263657 … 0.478193 0.448695 0.407768 0.421938 0.382411 0.355439 0.348323 0.532393 0.498619 0.46191
Display.
figure(figsize = (5,5))
imageplot(clamP(fM), @sprintf("Linear, Fourier, SNR = %.1f dB", snr(f, fM)))
PyObject <matplotlib.text.Text object at 0x0000000028CD61D0>
Given a $T$, the number of coefficients is obtained by counting the non-thresholded coefficients $ \abs{I_M} $.
m = sum(FT .!= 0)
print(@sprintf("M/N = 1/%d" ,(n0^2)/m))
M/N = 1/31
Exercise 2
Compute the value of the threshold $T$ so that the number of coefficients is $M$. Display the corresponding approximation $f_M$.
#run -i nt_solutions/introduction_4_fourier_wavelets/exo2
include("introduction_4_fourier_wavelets\ exo2.jl")
PyObject <matplotlib.text.Text object at 0x0000000028A56860>
## Insert your code here.
A wavelet basis $ \Bb = \{ \psi_m \}_m $ is obtained over the continuous domain by translating and dilating three mother wavelet functions $ \{\psi^V,\psi^H,\psi^D\} $.
Each wavelet atom is defined as $$ \psi_m(x) = \psi_{j,n}^k(x) = \frac{1}{2^j}\psi^k\pa{ \frac{x-2^j n}{2^j} } $$
The scale (size of the support) is $2^j$ and the position is $2^j(n_1,n_2)$. The index is $ m=(k,j,n) $ for $\{ j \leq 0 \}$.
The wavelet transform computes all the inner products $ \{ \dotp{f}{\psi_{j,n}^k} \}_{k,j,n} $.
Set the minimum scale for the transform to be 0.
Jmin = 0
0
Perform the wavelet transform, $f_w$ stores all the wavelet coefficients.
fw = NtToolBox.perform_wavelet_transf(f, Jmin, + 1)
#xf = f[1:2^(5+1), 1:2^(5+1)]
#vcat(xf[1:end, :], xf[end, :]')
#xf[1:end, :]
#xf[end, :]
512×512 Array{Float32,2}: 4.51785f10 1.90135f10 1.13727f9 … 4.80504 4.67925 3.88752 1.8792f10 7.91092f9 1.07963f9 4.69641 4.49987 3.67979 1.13615f9 1.17975f9 4.74278f8 4.23611 3.88736 3.05308 1.03742f9 1.10204f9 4.28281f8 2.88551 2.24264 1.65009 7.30674f7 7.25998f7 7.79994f7 1.44125 1.15025 0.999714 6.85523f7 6.7547f7 7.4514f7 … 0.984214 0.924601 0.910391 6.25727f7 5.92548f7 6.74838f7 0.786303 0.89966 0.9531 5.96832f7 5.55694f7 6.6863f7 0.826164 0.894299 0.950045 4.90871f6 4.51256f6 4.30903f6 0.823292 0.748641 0.951321 4.37187f6 4.21775f6 4.18829f6 0.808546 0.789063 1.0126 4.2607f6 4.16869f6 4.10386f6 … 0.85383 0.849685 1.0396 4.17215f6 4.0568f6 3.7211f6 0.865919 0.94616 1.16593 3.88938f6 3.84241f6 3.26856f6 0.982855 1.06794 1.31065 ⋮ ⋱ ⋮ 1.16031 1.61989 2.25564 … 0.772912 0.629928 0.534481 1.01383 1.31669 1.87845 0.655793 0.446958 0.366456 0.820397 1.02855 1.37966 0.498542 0.46958 0.391482 0.92509 1.09816 1.28413 0.481063 0.460606 0.467485 1.31634 1.29588 1.36797 0.429293 0.52516 0.511014 1.21311 1.24159 1.40131 … 0.410844 0.462497 0.545557 1.06359 1.11973 1.20602 0.53408 0.55677 0.634675 1.0181 1.08811 1.11144 0.642725 0.77148 0.94774 0.993192 1.09067 1.12538 0.825504 0.918117 1.01978 1.12733 1.038 1.08189 0.988706 1.0371 1.01898 0.85768 0.921585 0.994552 … 1.07007 1.1659 1.13605 0.910494 1.12769 1.18632 1.14715 1.2112 1.21211
using NPZ
test = npzread("pgsh.npy")
512×512 Array{Float64,2}: 233.032 19.1135 -13.0311 … 0.0418207 -0.131423 -23.1363 1.65645 -14.2758 0.0413636 -0.128621 -4.93937 15.4845 1.19021 0.0394703 -0.133461 -5.88592 -3.01506 -3.32998 -0.00103448 -0.0385822 0.922625 -2.16059 0.014336 -0.0159323 -0.00210343 -0.335541 5.64581 -9.66819 … -0.00869429 0.000518759 -0.518646 -0.502312 5.0826 0.00102269 0.0117643 4.98911 -5.50287 9.23286 -0.00281264 -0.00380574 0.964176 0.19277 0.00966503 -0.0331573 0.012915 -3.20033 0.128032 -0.512992 0.00189917 0.0302912 0.319147 0.377833 0.160483 … 0.0115657 0.0148096 -0.273935 -0.155208 0.166632 0.00603326 0.0116551 -0.603172 0.766374 -1.38124 -0.00246024 0.0250189 ⋮ ⋱ ⋮ -0.00091898 -0.0039071 -0.0059023 … -0.0110338 -0.00267938 0.0181912 -0.00505674 0.00412992 -0.0279641 -0.0170076 -0.0186527 -0.0144188 -0.0191199 0.00409298 -0.00667346 -0.0275172 0.00556199 -0.0102931 -0.0184985 0.000831558 0.0513395 -0.00900993 0.00858931 0.018305 0.0116154 0.0127468 -0.0135157 0.027699 … -0.00627344 0.0268372 0.00047219 -0.0162774 -0.0247788 -0.00592464 -0.0204156 0.00738821 -4.32492e-5 -0.00901067 0.011207 0.0168448 -0.0102837 0.00511243 -0.00709129 -0.00368487 -0.0161069 0.0455343 -0.0165399 -0.000882116 0.00323547 -0.0198884 -0.010534 -0.0135025 -0.00207893 … 0.00426511 -0.00434588 2.50898e-5 0.0201377 0.0201784 -3.24896e-5 -0.000651598
Display the transformed coefficients.
figure(figsize = (10, 10))
NtToolBox.plot_wavelet(fw)
title("Wavelet coefficients")
PyObject <matplotlib.text.Text object at 0x0000000028837278>
Linear wavelet approximation with $M=2^{-j_0}$ coefficients is obtained by keeping only the coarse scale (large support) wavelets:
$$ I_M = \enscond{(k,j,n)}{ j \geq j_0 }. $$It corresponds to setting to zero all the coefficients excepted those that are on the upper left corner of $f_w$.
Exercise 3
Perform linear approximation with $M$ wavelet coefficients.
run -i nt_solutions/introduction_4_fourier_wavelets/exo3
## Insert your code here.
A non-linear approximation is obtained by keeping the $M$ largest wavelet coefficients.
As already said, this is equivalently computed by a non-linear hard thresholding.
Select a threshold.
T = .15
0.15
Perform hard thresholding.
fwT = fw*(abs(fw) .> T)
512×512 Array{Float32,2}: 6.68709f10 6.68709f10 6.68709f10 … 6.68709f10 6.68709f10 2.93133f10 2.93133f10 2.93133f10 2.93133f10 2.93133f10 3.61903f9 3.61903f9 3.61903f9 3.61903f9 3.61903f9 3.37039f9 3.37039f9 3.37039f9 3.37039f9 3.37039f9 4.60282f8 4.60282f8 4.60282f8 4.60282f8 4.60282f8 4.47582f8 4.47582f8 4.47582f8 … 4.47582f8 4.47582f8 4.16612f8 4.16612f8 4.16612f8 4.16612f8 4.16612f8 4.04895f8 4.04895f8 4.04895f8 4.04895f8 4.04895f8 5.71657f7 5.71657f7 5.71657f7 5.71657f7 5.71657f7 5.65475f7 5.65475f7 5.65475f7 5.65475f7 5.65475f7 5.68546f7 5.68546f7 5.68546f7 … 5.68546f7 5.68546f7 5.22421f7 5.22421f7 5.22421f7 5.22421f7 5.22421f7 5.00265f7 5.00265f7 5.00265f7 5.00265f7 5.00265f7 ⋮ ⋱ ⋮ 1110.42 1110.42 1110.42 … 1109.89 1110.42 1112.5 1112.5 1112.5 1112.06 1112.5 1110.06 1110.06 1110.06 1109.63 1110.06 1122.24 1122.24 1122.24 1121.76 1122.24 1127.55 1127.55 1127.55 1127.01 1127.55 1129.02 1129.02 1129.02 … 1128.43 1129.02 1139.33 1139.33 1139.33 1138.85 1139.33 1160.46 1160.46 1160.46 1160.12 1160.46 1167.48 1167.48 1167.48 1167.12 1167.48 1174.75 1174.75 1174.75 1174.3 1174.75 1176.6 1176.6 1176.6 … 1176.13 1176.6 1179.26 1179.26 1179.26 1178.83 1179.26
Display the thresholded coefficients.
figure(figsize = (15,15))
subplot(1, 2, 1)
plot_wavelet(fw)
title("Original coefficients")
subplot(1, 2, 2)
plot_wavelet(fwT)
title("Thresholded coefficients")
show()
Perform reconstruction.
fM = perform_wavelet_transf(fwT, Jmin, -1)
# x = copy(fwT)
# a = transpose(x[:, 1])
# fin = length(a)
# d = x[Int(fin/2) + 1: end, :].*h[end]
# x = x[1:Int(fin/2), :]./h[end]
ArgumentError: invalid index: 2.0 in macro expansion at .\multidimensional.jl:350 [inlined] in macro expansion at .\cartesian.jl:64 [inlined] in macro expansion at .\multidimensional.jl:348 [inlined] in _unsafe_getindex! at .\multidimensional.jl:340 [inlined] in macro expansion at .\multidimensional.jl:298 [inlined] in _unsafe_getindex(::Base.LinearFast, ::Array{Float32,2}, ::FloatRange{Float64}, ::Colon) at .\multidimensional.jl:291 in lifting_step(::Array{Float32,2}, ::Array{Float64,1}, ::Int64) at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:178 in perform_wavelet_transf(::Array{Float32,2}, ::Int64, ::Int64, ::String, ::Int64, ::Int64) at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:90 in perform_wavelet_transf(::Array{Float32,2}, ::Int64, ::Int64) at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:3
Display approximation.
plt.figure(figsize=(5,5))
imageplot(clamp(fM), "Approximation, SNR, = %.1f dB" %snr(f, fM))
Exercise 4
Perform non-linear approximation with $M$ wavelet coefficients by chosing the correct value for $T$. Store the result in the variable $f_M$.
run -i nt_solutions/introduction_4_fourier_wavelets/exo4
## Insert your code here.
Compare two 1D profile (lines of the image). Note how the ringing artifacts are reduced compared to the Fourier approximation.
plt.figure(figsize=(7,6))
plt.subplot(2, 1, 1)
plt.plot(f[:,n0//2])
plt.title('f')
plt.subplot(2, 1, 2)
plt.plot(fM[:,n0//2])
plt.title('f_M')
plt.show()