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 explores 2-D multiresolution analysis with Daubechies wavelet transform.
# using Autoreload
# arequire("nt_toolbox/nt_signal.jl")
# arequire("nt_toolbox/nt_general.jl")
# using nt_signal
# using nt_general
using PyPlot
using NtToolBox
The 2-D wavelet transform of a continuous image $f(x)$ computes the set of inner products $$ d_j^k[n] = \dotp{f}{\psi_{j,n}^k} $$ for scales $ j \in \ZZ $, position $ n \in \ZZ^2 $ and orientation $ k \in \{H,V,D\} $.
The wavelet atoms are defined by scaling and translating three mother atoms $ \{\psi^H,\psi^V,\psi^D\} $: $$ \psi_{j,n}^k(x) = \frac{1}{2^j}\psi^k\pa{\frac{x-2^j n}{2^j}} $$ These oriented wavelets are defined by a tensor product of a 1-D wavelet function $\psi(t)$ and a 1-D scaling function $\phi(t)$ $$ \psi^H(x)=\phi(x_1)\psi(x_2), \quad \psi^V(x)=\psi(x_1)\phi(x_2) \qandq \psi^D(x)=\psi(x_1)\psi(x_2).$$
The fast wavelet transform algorithm does not make use of the wavelet and scaling functions, but of the filters $h$ and $g$ that caracterize their interaction: $$ g[n] = \frac{1}{\sqrt{2}}\dotp{\psi(t/2)}{\phi(t-n)} \qandq h[n] = \frac{1}{\sqrt{2}}\dotp{\phi(t/2)}{\phi(t-n)}. $$
The simplest filters are the Haar filters $$ h = [1, 1]/\sqrt{2} \qandq g = [-1, 1]/\sqrt{2}. $$
Daubechies wavelets extends the haar wavelets by using longer filters, that produce smoother scaling functions and wavelets. Furthermore, the larger the size $p=2k$ of the filter, the higher is the number $k$ of vanishing moment.
A high number of vanishing moments allows to better compress regular parts of the signal. However, increasing the number of vanishing moments also inceases the size of the support of the wavelets, wich can be problematic in part where the signal is singular (for instance discontinuous).
Choosing the best wavelet, and thus choosing $k$, that is adapted to a given class of signals, thus corresponds to a tradeoff between efficiency in regular and singular parts.
Set the support size. To begin, we select the D4 filter.
p = 4;
Create the low pass filter $h$ and the high pass $g$. We add a zero to ensure that it has a odd length. Note that the central value of $h$ corresponds to the 0 position.
h = [0, .482962913145, .836516303738, .224143868042, -.129409522551];
h = h/norm(h);
g = cat(1, 0, h[length(h):-1:2]) .* ( (-1).^(1:length(h)) );
Note that the high pass filter $g$ is computed directly from the low pass filter as: $$g[n] = (-1)^{1-n}h[1-n]$$
Display.
println("h filter = ", h);
println("g filter = ", g);
h filter = [0.0,0.482963,0.836516,0.224144,-0.12941] g filter = [-0.0,-0.12941,-0.224144,0.836516,-0.482963]
The basic wavelet operation is low/high filtering, followed by down sampling.
Starting from some 1-D signal $f \in \RR^N$, one thus compute the low pass signal $a \in \RR^{N/2}$ and the high pass signal $d \in \RR^{N/2}$ as $$ a = (f \star h) \downarrow 2 \qandq d = (f \star g) \downarrow 2$$ where the sub-sampling is defined as $$ (u \downarrow 2)[k] = u[2k]. $$
Create a random signal $f \in \RR^N$.
N = 256;
f = rand(N,1);
Low/High pass filtering followed by sub-sampling.
a = subsampling( cconvol(f,h) );
d = subsampling( cconvol(f,g) );
UndefVarError: subsampling not defined
For orthogonal filters, the reverse of this process is its dual (aka its transpose), which is upsampling followed by low/high pass filtering with the reversed filters and summing: $$ (a \uparrow h) \star \tilde h + (d \uparrow g) \star \tilde g = f $$ where $\tilde h[n]=h[-n]$ (computed modulo $N$) and $ (u \uparrow 2)[2n]=u[n] $ and $ (u \uparrow 2)[2n+1]=0 $.
Test for energy conservation.
e0 = norm(f[:])^2;
e1 = norm(a[:])^2 + norm(d[:])^2;
println("Energy before : $e0");
println("Energy before : $e1");
UndefVarError: a not defined
Up-sampling followed by filtering.
f1 = cconvol(upsampling(a),reverse(h)) + cconvol(upsampling(d),reverse(g));
UndefVarError: cconvol not defined
Check that we really recover the same signal.
println("Error |f-f1|/|f| = ", norm(f[:]-f1[:])/norm(f[:]) );
UndefVarError: f1 not defined
The set of wavelet coefficients are computed with a fast algorithm that exploit the embedding of the approximation spaces $V_j$ spanned by the scaling function $ \{ \phi_{j,n} \}_n $ defined as $$ \phi_{j,n}(x) = \frac{1}{2^j}\phi^0\pa{\frac{x-2^j n}{2^j}} \qwhereq \phi^0(x)=\phi(x_1)\phi(x_2). $$
The wavelet transform of $f$ is computed by using intermediate discretized low resolution images obtained by projection on the spaces $V_j$: $$ a_j[n] = \dotp{f}{\phi_{j,n}}. $$
First we load an image of $N= n \times n$ pixels.
n = 256;
name = "nt_toolbox/data/hibiscus.bmp";
f = load_image(name, N);
f = rescale(sum(f,3));
f = f[:,:,1];
imageplot(f);
PyError (:PyObject_Call) <type 'exceptions.ValueError'> ValueError(u"Only know how to handle extensions: [u'png']; with Pillow installed matplotlib can handle more images",) File "/Users/gpeyre/.julia/v0.5/Conda/deps/usr/lib/python2.7/site-packages/matplotlib/pyplot.py", line 2314, in imread return _imread(*args, **kwargs) File "/Users/gpeyre/.julia/v0.5/Conda/deps/usr/lib/python2.7/site-packages/matplotlib/image.py", line 1282, in imread 'more images' % list(six.iterkeys(handlers))) in pyerr_check at /Users/gpeyre/.julia/v0.5/PyCall/src/exception.jl:56 [inlined] in pyerr_check at /Users/gpeyre/.julia/v0.5/PyCall/src/exception.jl:61 [inlined] in macro expansion at /Users/gpeyre/.julia/v0.5/PyCall/src/exception.jl:81 [inlined] in #_pycall#66(::Array{Any,1}, ::Function, ::PyCall.PyObject, ::String, ::Vararg{String,N}) at /Users/gpeyre/.julia/v0.5/PyCall/src/PyCall.jl:620 in _pycall(::PyCall.PyObject, ::String, ::Vararg{String,N}) at /Users/gpeyre/.julia/v0.5/PyCall/src/PyCall.jl:608 in #pycall#70(::Array{Any,1}, ::Function, ::PyCall.PyObject, ::Type{PyCall.PyAny}, ::String, ::Vararg{String,N}) at /Users/gpeyre/.julia/v0.5/PyCall/src/PyCall.jl:642 in pycall(::PyCall.PyObject, ::Type{PyCall.PyAny}, ::String, ::Vararg{String,N}) at /Users/gpeyre/.julia/v0.5/PyCall/src/PyCall.jl:642 in #imread#65(::Array{Any,1}, ::Function, ::String, ::Vararg{String,N}) at /Users/gpeyre/.julia/v0.5/PyPlot/src/PyPlot.jl:172 in load_image(::String, ::Int64, ::Int64, ::Int64, ::Int64) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:12 (repeats 2 times)
The algorithm starts at the coarsest scale $ j=\log_2(n)-1 $
j = log2(n)-1;
The first step of the algorithm perform filtering/downsampling in the horizontal direction.
$$ \tilde a_{j-1} = (a_j \star^H h) \downarrow^{2,H} \qandq \tilde d_{j-1} = (a_j \star^H g) \downarrow^{2,H}$$Here, the operator $\star^H$ and $\downarrow^{2,H}$ are defined by applying $\star$ and $\downarrow^2$ to each column of the matrix.
The second step computes the filtering/downsampling in the vertical direction.
$$ a_{j-1} = (\tilde a_j \star^V h) \downarrow^{2,V} \qandq d_{j-1}^V = (\tilde a_j \star^V g) \downarrow^{2,V},$$$$ d_{j-1}^H = (\tilde d_j \star^V h) \downarrow^{2,V} \qandq d_{j-1}^D = (\tilde d_j \star^V g) \downarrow^{2,V}.$$A wavelet transform is computed by iterating high pass and loss pass filterings with |h| and |g|, followed by sub-samplings. Since we are in 2-D, we need to compute these filterings+subsamplings in the horizontal and then in the vertical direction (or in the reverse order, it does not mind).
Initialize the transformed coefficients as the image itself and set the initial scale as the maximum one. |fW| will be iteratively transformated and will contains the coefficients.
fW = copy(f);
Select the sub-part of the image to transform.
A = fW[1:2^(j+1),1:2^(j+1)];
BoundsError: attempt to access 256×1 Array{Float64,2} at index [1.0:1.0:256.0,1.0:1.0:256.0] in throw_boundserror(::Array{Float64,2}, ::Tuple{FloatRange{Float64},FloatRange{Float64}}) at ./abstractarray.jl:355 in checkbounds at ./abstractarray.jl:284 [inlined] in _getindex at ./multidimensional.jl:270 [inlined] in getindex(::Array{Float64,2}, ::FloatRange{Float64}, ::FloatRange{Float64}) at ./abstractarray.jl:752
Apply high and low filtering+subsampling in the vertical direction (1st ooordinate), to get coarse and details.
Coarse = subsampling(cconvol(A,h,1),1);
Detail = subsampling(cconvol(A,g,1),1);
UndefVarError: subsampling not defined
Check for energy conservation.
norm(A[:])^2 - norm(Coarse[:])^2 - norm(Detail[:])^2
UndefVarError: A not defined
Note: |subsamplling(A,1)| is equivalent to |A(1:2:end,:)| and |subsamplling(A,2)| is equivalent to |A(:,1:2:end)|.
Concatenate them in the vertical direction to get the result.
A = cat(1, Coarse, Detail );
UndefVarError: Coarse not defined
Display the result of the vertical transform.
clf;
imageplot(f,"Original image",1,2,1);
imageplot(A,"Vertical transform",1,2,2);
MethodError: no method matching imageplot(::Array{Float64,2}, ::String, ::Int64, ::Int64, ::Int64) Closest candidates are: imageplot(::Any, ::Any, ::Any) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:45 imageplot(::Any, ::Any) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:45 imageplot(::Any) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:45
Apply high and low filtering+subsampling in the horizontal direction (2nd ooordinate), to get coarse and details.
Coarse = subsampling(cconvol(A,h,2),2);
Detail = subsampling(cconvol(A,g,2),2);
UndefVarError: subsampling not defined
Concatenate them in the horizontal direction to get the result.
A = cat(2, Coarse, Detail );
UndefVarError: Coarse not defined
Assign the transformed data.
fW[1:2^(j+1),1:2^(j+1)] = A;
UndefVarError: A not defined
Display the result of the horizontal transform.
clf;
imageplot(f,"Original image",1,2,1);
subplot(1,2,2);
plot_wavelet(fW,log2(n)-1);
title("Transformed");
MethodError: no method matching imageplot(::Array{Float64,2}, ::String, ::Int64, ::Int64, ::Int64) Closest candidates are: imageplot(::Any, ::Any, ::Any) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:45 imageplot(::Any, ::Any) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:45 imageplot(::Any) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:45
print( norm(f[:])-norm(fW[:]) )
0.0
Exercise 1
Implement a full wavelet transform that extract iteratively wavelet coefficients, by repeating these steps. Take care of choosing the correct number of steps.
include("nt_solutions/wavelet_4_daubechies2d/exo1.jl");
LoadError: BoundsError: attempt to access 256×1 Array{Float64,2} at index [1.0:1.0:256.0,1.0:1.0:256.0] while loading /Users/gpeyre/Dropbox/github/numerical-tours/julia/nt_solutions/wavelet_4_daubechies2d/exo1.jl, in expression starting on line 5 in throw_boundserror(::Array{Float64,2}, ::Tuple{FloatRange{Float64},FloatRange{Float64}}) at ./abstractarray.jl:355 in checkbounds at ./abstractarray.jl:284 [inlined] in _getindex at ./multidimensional.jl:270 [inlined] in getindex(::Array{Float64,2}, ::FloatRange{Float64}, ::FloatRange{Float64}) at ./abstractarray.jl:752 in macro expansion; at /Users/gpeyre/Dropbox/github/numerical-tours/julia/nt_solutions/wavelet_4_daubechies2d/exo1.jl:6 [inlined] in anonymous at ./<missing>:? in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
Check for orthogonality of the transform (conservation of energy).
e0 = norm(f[:])^2;
e1 = norm(fW[:])^2;
println("Energy of the signal = $e0");
println("Energy of the coefficients = $e1");
Energy of the signal = 81.92047340519909 Energy of the coefficients = 81.92047340519909
Display the wavelet coefficients.
clf;
imageplot(f, "Original", 1,2,1);
subplot(1,2,2);
plot_wavelet(fW, Jmin);
title("Transformed");
MethodError: no method matching imageplot(::Array{Float64,2}, ::String, ::Int64, ::Int64, ::Int64) Closest candidates are: imageplot(::Any, ::Any, ::Any) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:45 imageplot(::Any, ::Any) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:45 imageplot(::Any) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:45
Inversing the wavelet transform means retrieving a signal |f1| from the coefficients |fW|. If |fW| are exactely the coefficients of |f|, then |f=f1| up to machine precision.
Initialize the image to recover |f1| as the transformed coefficient, and select the smallest possible scale.
f1 = copy(fW);
j = 0;
Select the sub-coefficient to transform.
A = f1[1:2^(j+1),1:2^(j+1)];
BoundsError: attempt to access 256×1 Array{Float64,2} at index [1:2,1:2] in throw_boundserror(::Array{Float64,2}, ::Tuple{UnitRange{Int64},UnitRange{Int64}}) at ./abstractarray.jl:355 in checkbounds at ./abstractarray.jl:284 [inlined] in _getindex at ./multidimensional.jl:270 [inlined] in getindex(::Array{Float64,2}, ::UnitRange{Int64}, ::UnitRange{Int64}) at ./abstractarray.jl:752
Retrieve coarse and detail coefficients in the vertical direction (you can begin by the other direction, this has no importance).
Coarse = A[1:2^j,:];
Detail = A[2^j+1:2^(j+1),:];
UndefVarError: A not defined
Undo the transform by up-sampling and then dual filtering.
Coarse = cconvol(upsampling(Coarse,1),reverse(h),1);
Detail = cconvol(upsampling(Detail,1),reverse(g),1);
Recover the coefficient by summing.
A = Coarse + Detail;
UndefVarError: Coarse not defined
Retrieve coarse and detail coefficients in the vertical direction (you can begin by the other direction, this has no importance).
Coarse = A[:,1:2^j];
Detail = A[:,2^j+1:2^(j+1)];
UndefVarError: A not defined
Undo the transform by up-sampling and then dual filtering.
Coarse = cconvol(upsampling(Coarse,2),reverse(h),2);
Detail = cconvol(upsampling(Detail,2),reverse(g),2);
UndefVarError: cconvol not defined
Recover the coefficient by summing.
A = Coarse + Detail;
UndefVarError: Coarse not defined
Assign the result.
f1[1:2^(j+1),1:2^(j+1)] = A;
UndefVarError: A not defined
Exercise 2
Write the inverse wavelet transform that computes |f1| from the coefficients |fW|. Compare |f1| with |f|.
include("nt_solutions/wavelet_4_daubechies2d/exo2.jl");
LoadError: BoundsError: attempt to access 256×1 Array{Float64,2} at index [1.0:1.0:4.0,1.0:1.0:4.0] while loading /Users/gpeyre/Dropbox/github/numerical-tours/julia/nt_solutions/wavelet_4_daubechies2d/exo2.jl, in expression starting on line 3 in throw_boundserror(::Array{Float64,2}, ::Tuple{FloatRange{Float64},FloatRange{Float64}}) at ./abstractarray.jl:355 in checkbounds at ./abstractarray.jl:284 [inlined] in _getindex at ./multidimensional.jl:270 [inlined] in getindex(::Array{Float64,2}, ::FloatRange{Float64}, ::FloatRange{Float64}) at ./abstractarray.jl:752 in macro expansion; at /Users/gpeyre/Dropbox/github/numerical-tours/julia/nt_solutions/wavelet_4_daubechies2d/exo2.jl:4 [inlined] in anonymous at ./<missing>:? in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
Check that we recover exactly the original image.
e = norm(f[:]-f1[:])/norm(f[:]);
print("Error |f-f1|/|f| = $e");
Error |f-f1|/|f| = 0.0
Linear approximation is performed by setting to zero the fine scale wawelets coefficients and then performing the inverse wavelet transform.
Here we keep only 1/16 of the wavelet coefficient, thus calculating an |m| term approximation with |m=n^2/16|.
eta = 4;
fWLin = zeros(n,n);
fWLin[1:n/eta,1:n/eta] = fW[1:n/eta,1:n/eta];
BoundsError: attempt to access 256×1 Array{Float64,2} at index [1.0:1.0:64.0,1.0:1.0:64.0] in throw_boundserror(::Array{Float64,2}, ::Tuple{FloatRange{Float64},FloatRange{Float64}}) at ./abstractarray.jl:355 in checkbounds at ./abstractarray.jl:284 [inlined] in _getindex at ./multidimensional.jl:270 [inlined] in getindex(::Array{Float64,2}, ::FloatRange{Float64}, ::FloatRange{Float64}) at ./abstractarray.jl:752
Exercise 3: Compute and display the linear approximation |fLin| obtained from the coefficients |fWLin| by inverse wavelet transform.
include("nt_solutions/wavelet_4_daubechies2d/exo3.jl");
LoadError: BoundsError: attempt to access 256×1 Array{Float64,2} at index [1:256,1:256] while loading /Users/gpeyre/Dropbox/github/numerical-tours/julia/nt_solutions/wavelet_4_daubechies2d/exo3.jl, in expression starting on line 3 in throw_boundserror(::Array{Float64,2}, ::Tuple{UnitRange{Int64},UnitRange{Int64}}) at ./abstractarray.jl:355 in checkbounds at ./abstractarray.jl:284 [inlined] in _getindex at ./multidimensional.jl:270 [inlined] in getindex at ./abstractarray.jl:752 [inlined] in subselect(::Array{Float64,2}, ::UnitRange{Int64}) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:333 in perform_wavortho_transf(::Array{Float64,2}, ::Int64, ::Int64, ::Array{Float64,1}) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:364 in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?
A non-linear |m|-term approximation is obtained by keeping only the |m| largest coefficients, which creates the smallest possible error.
Removing the smallest coefficient, to keep the |m|-largest, is equivalently obtainedby thresholding the coefficients to set to 0 the smallest coefficients.
First select a threshold value (the largest the threshold, the more agressive the approximation).
T = .2;
Then set to 0 coefficients with magnitude below the threshold.
fWT = fW .* float(abs(fW).>T);
Display thresholded coefficients.
clf;
subplot(1,2,1);
plot_wavelet(fW,Jmin);
title("Original coefficients");
subplot(1,2,2);
plot_wavelet(fWT,Jmin);
title("Thresholded coefficients");
BoundsError: attempt to access 256×1 Array{Float64,2} at index [1:2,1:2] in throw_boundserror(::Array{Float64,2}, ::Tuple{UnitRange{Int64},UnitRange{Int64}}) at ./abstractarray.jl:355 in checkbounds at ./abstractarray.jl:284 [inlined] in _getindex at ./multidimensional.jl:270 [inlined] in getindex at ./abstractarray.jl:752 [inlined] in plot_wavelet(::Array{Float64,2}, ::Int64) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:278
Exercise 4
Find the thresholds |T| so that the number |m| of remaining coefficients in |fWT| are |m=n^2/16|. Use this threshold to compute |fWT| and then display the corresponding approximation |Mnlin| of |f|. Compare this result with the linear approximation. umber of kept coefficients ompute the threshold T elect threshold nverse transform nverse isplay
include("nt_solutions/wavelet_4_daubechies2d/exo4.jl");
LoadError: BoundsError: attempt to access 256×1 Array{Float64,2} at index [1:256,1:256] while loading /Users/gpeyre/Dropbox/github/numerical-tours/julia/nt_solutions/wavelet_4_daubechies2d/exo4.jl, in expression starting on line 5 in throw_boundserror(::Array{Float64,2}, ::Tuple{UnitRange{Int64},UnitRange{Int64}}) at ./abstractarray.jl:355 in checkbounds at ./abstractarray.jl:284 [inlined] in _getindex at ./multidimensional.jl:270 [inlined] in getindex at ./abstractarray.jl:752 [inlined] in subselect(::Array{Float64,2}, ::UnitRange{Int64}) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:333 in perform_wavortho_transf(::Array{Float64,2}, ::Int64, ::Int64, ::Array{Float64,1}) at /Users/gpeyre/.julia/v0.5/NtToolBox/src/signal.jl:364 in include_from_node1(::String) at ./loading.jl:488 in include_from_node1(::String) at /Applications/Julia-0.5.app/Contents/Resources/julia/lib/julia/sys.dylib:?