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 presents block thresholding methods, that makes use of the structure of wavelet coefficients of natural images to perform denoising. Theoretical properties of block thresholding were investigated in CaiSilv Cai99 HallKerkPic99
using PyPlot
using NtToolBox
using Autoreload
arequire("NtToolBox")
WARNING: Base.UTF8String is deprecated, use String instead.
likely near C:\Users\Ayman\.julia\v0.5\Autoreload\src\constants.jl:1
WARNING: Base.UTF8String is deprecated, use String instead.
likely near C:\Users\Ayman\.julia\v0.5\Autoreload\src\constants.jl:1
WARNING: Base.UTF8String is deprecated, use String instead.
likely near C:\Users\Ayman\.julia\v0.5\Autoreload\src\constants.jl:7
WARNING: Base.UTF8String is deprecated, use String instead.
likely near C:\Users\Ayman\.julia\v0.5\Autoreload\src\constants.jl:7
WARNING: Base.UTF8String is deprecated, use String instead.
likely near C:\Users\Ayman\.julia\v0.5\Autoreload\src\constants.jl:12
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: Base.UTF8String is deprecated, use String instead.
likely near In[1]:4
WARNING: readall is deprecated, use readstring instead.
in depwarn(::String, ::Symbol) at .\deprecated.jl:64
in readall(::IOStream, ::Vararg{IOStream,N}) at .\deprecated.jl:30
in #parse_file#6(::Array{Any,1}, ::Function, ::String) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\files.jl:63
in parse_file(::String) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\files.jl:61
in #arequire#10(::Symbol, ::Array{String,1}, ::Function, ::String) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:65
in arequire(::String) at C:\Users\Ayman\.julia\v0.5\Autoreload\src\Autoreload.jl:47
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[1], in expression starting on line 4
WARNING: replacing module NtToolBox
Here we use an additive Gaussian noise.
Size of the image of $N=n \times n$ pixels.
n = 256
WARNING: Base.UTF8String is deprecated, use String instead. likely near C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl:31 WARNING: Base.UTF8String is deprecated, use
256
String instead. likely near C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl:31 WARNING: Base.UTF8String is deprecated, use String instead. likely near C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl:31 WARNING: Base.UTF8String is deprecated, use String instead. likely near C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl:31 WARNING: Base.UTF8String is deprecated, use String instead. likely near C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl:31 WARNING: Base.UTF8String is deprecated, use String instead. likely near C:\Users\Ayman\.julia\v0.5\IJulia\src\kernel.jl:31
First we load an image $f_0 \in \RR^N$.
f0 = rescale(load_image("NtToolBox/src/data/boat.png", n))
256×256 Array{Float32,2}: 0.726186 0.726186 0.730672 0.726341 … 0.585907 0.590307 0.572848 0.726169 0.730517 0.726306 0.721822 0.590402 0.590298 0.603375 0.713039 0.721684 0.726118 0.726135 0.586037 0.581663 0.585985 0.708606 0.69962 0.708605 0.712937 0.590471 0.59469 0.58608 0.695716 0.68678 0.699757 0.708536 0.599319 0.616285 0.581723 0.699721 0.691299 0.695833 0.699756 … 0.611987 0.594681 0.590693 0.712834 0.704172 0.682702 0.691332 0.607716 0.60771 0.599192 0.704155 0.695765 0.699555 0.703816 0.595651 0.620214 0.611465 0.674147 0.69099 0.695139 0.691042 0.603527 0.603503 0.603608 0.686745 0.678441 0.673822 0.695308 0.616428 0.611733 0.599415 0.687088 0.686932 0.674483 0.678379 … 0.604203 0.604266 0.612237 0.673787 0.678408 0.69488 0.686783 0.620736 0.608477 0.612091 0.682774 0.68671 0.694801 0.686292 0.624694 0.620577 0.617054 ⋮ ⋱ ⋮ 0.645094 0.431234 0.354052 0.308235 0.494348 0.53709 0.403382 0.690225 0.600543 0.636972 0.744182 … 0.494731 0.405744 0.374266 0.608841 0.608419 0.629718 0.601788 0.459566 0.431199 0.46154 0.675206 0.743989 0.779907 0.767567 0.336347 0.521763 0.479617 0.714329 0.667938 0.633161 0.6641 0.440696 0.380352 0.384329 0.612804 0.659272 0.677521 0.680758 0.594358 0.50617 0.531426 0.750684 0.71642 0.647809 0.661189 … 0.620309 0.494524 0.424604 0.700667 0.717274 0.681865 0.773526 0.478253 0.580802 0.389469 0.690826 0.673432 0.565668 0.65705 0.515763 0.542494 0.475566 0.673987 0.659925 0.65977 0.769325 0.495179 0.446196 0.515317 0.603688 0.651824 0.74315 0.695149 0.497986 0.612426 0.267147 0.682142 0.682239 0.638636 0.656232 … 0.383471 0.481013 0.0
Display it.
figure(figsize = (5,5))
imageplot(f0)
Noise level.
sigma = .08
0.08
Generate a noisy image $f=f_0+\epsilon$ where $\epsilon \sim \Nn(0,\si^2\text{Id}_N)$.
using Distributions
f = f0 + sigma.*rand(Normal(), n, n)
256×256 Array{Float64,2}: 0.710313 0.793932 0.852947 0.712922 … 0.64466 0.463513 0.560563 0.710355 0.715417 0.692543 0.796806 0.584508 0.720498 0.652994 0.768987 0.685349 0.564781 0.817171 0.649122 0.589976 0.640439 0.656056 0.691161 0.609962 0.730005 0.694955 0.58525 0.665006 0.739125 0.731892 0.891632 0.701638 0.50126 0.667541 0.520729 0.803983 0.923905 0.783055 0.613976 … 0.661443 0.519405 0.596635 0.636989 0.819037 0.672673 0.733444 0.591171 0.614601 0.542568 0.693603 0.813594 0.733113 0.647362 0.549997 0.679936 0.572804 0.662841 0.738771 0.653678 0.752576 0.468689 0.527867 0.736143 0.860663 0.6223 0.683235 0.729511 0.582323 0.71925 0.545089 0.711831 0.700879 0.712673 0.716136 … 0.60234 0.612328 0.663829 0.705548 0.588976 0.694865 0.623223 0.619137 0.736872 0.615701 0.708127 0.552672 0.682251 0.565554 0.679682 0.538843 0.588373 ⋮ ⋱ ⋮ 0.606097 0.354045 0.341769 0.251088 0.49171 0.515805 0.392602 0.694127 0.720858 0.552058 0.686459 … 0.511795 0.357301 0.354399 0.529064 0.634521 0.688368 0.530925 0.358881 0.316357 0.428326 0.737423 0.765996 0.862529 0.802569 0.406906 0.520878 0.41159 0.707312 0.500645 0.733974 0.617463 0.472911 0.241862 0.224887 0.723371 0.653345 0.678984 0.744956 0.625732 0.376655 0.739158 0.801612 0.677163 0.543218 0.602471 … 0.545878 0.521769 0.390054 0.736641 0.776965 0.682361 0.747757 0.458962 0.615079 0.353581 0.614217 0.737972 0.575287 0.602417 0.560156 0.625889 0.43057 0.670669 0.543122 0.608007 0.718555 0.581604 0.453024 0.61008 0.72397 0.674379 0.910689 0.585294 0.496575 0.646263 0.215787 0.699264 0.702948 0.749954 0.51241 … 0.536127 0.390011 0.149293
Display it.
figure(figsize = (5,5))
imageplot(clamP(f))
We first consider the traditional wavelet thresholding method.
Parameters for the orthogonal wavelet transform.
Jmin = 4
4
Shortcuts for the foward and backward wavelet transforms.
wav = f -> perform_wavelet_transf(f, Jmin, +1)
iwav = fw -> perform_wavelet_transf(fw, Jmin, -1)
(::#3) (generic function with 1 method)
Display the original set of noisy coefficients.
figure(figsize = (10, 10))
plot_wavelet(wav(f), Jmin)
show()
Denoting $\Ww$ and $\Ww^*$ the forward and backward wavelet transform, wavelet thresholding $\tilde f$ is defined as
$$ \tilde f = \Ww^* \circ \theta_T \circ \Ww(f) $$where $T>0$ is the threshold, that should be adapted to the noise level.
The thresholding operator is applied component-wise
$$ \th_T(x)_i = \psi_T(x_i) x_i $$where $\psi_T$ is an atenuation fonction. In this tour, we use the James Stein (JS) attenuation:
$$ \psi_T(s) = \max\pa{ 0, 1-\frac{T^2}{s^2} } $$# psi = (s, T) -> max(1 - T^2./max(abs(s).^2, 1e-9.*ones(size(s)[1])), zeros(size(s)[1]))
psi = (s, T) -> max(1 - T^2./max(abs(s).^2, 1e-9.*ones(size(s))), zeros(size(s)))
(::#101) (generic function with 1 method)
Display the thresholding function $\th_T$.
s = linspace(-3, 3, 1024)
plot(s, s.*psi(s, 1))
plot(s, s, "r--")
show()
Thresholding operator.
theta = (x, T) -> psi(x, T).*x
ThreshWav = (f, T) -> iwav(theta(wav(f), T))
(::#15) (generic function with 1 method)
Test the thresholding.
T = 1.5*sigma
figure(figsize = (5, 5))
imageplot(clamP(ThreshWav(f, T)))
# I must find a solution for perform_wavelet_transf, the positive portion works but not well while the negative portion does not work
DimensionMismatch("dimensions must match") in promote_shape(::Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}, ::Tuple{Base.OneTo{Int64}}) at .\operators.jl:411 in indices(::Base.Zip2{Array{Float64,2},Array{Float64,1}}) at .\iterator.jl:94 in _array_for(::Type{Float64}, ::Base.Zip2{Array{Float64,2},Array{Float64,1}}, ::Base.HasShape) at .\array.jl:295 in collect(::Base.Generator{Base.Zip2{Array{Float64,2},Array{Float64,1}},Base.Math.##27#30}) at .\array.jl:308 in (::##11#12)(::Array{Float64,2}, ::Float64) at .\In[34]:1 in (::##13#14)(::Array{Float64,2}, ::Float64) at .\In[36]:1 in (::##15#16)(::Array{Float64,2}, ::Float64) at .\In[36]:2
Exercise 1
Display the evolution of the denoising SNR when $T$ varies. Store in $f_{Thresh}$ the optimal denoising result.
run -i nt_solutions/denoisingwav_4_block/exo1
## Insert your code here.
Display the optimal thresolding.
plt.figure(figsize = (5,5))
imageplot(clamp(fThresh), "SNR = %.1f dB" %snr(f0, fThresh))
A block thresholding operator of coefficients $x=(x_i)_{i=1}^P \in \RR^P$ is defined using a partition $B$ into a set of blocks $b$
$$ \{1,\ldots,P\} = \bigcup_{b \in B} b. $$Its definition reads
$$ \forall i \in b, \quad \theta_T(x)_i = \psi_T\left( \norm{x_b}_2 \right) x_i $$where $ x_b = (x_j)_{j \in B} \in \RR^{\abs{b}} $. One thus thresholds the $\ell^2$ norm (the energy) of each block rather than each coefficient independently.
For image-based thresholding, we use a partition in square blocks of equal size $w \times w$.
The block size $w$.
w = 4
n = 256
256
Compute indexing of the blocks.
include("ndgrid.jl")
(X, Y, dX, dY) = ndgrid(1:w:n-w+1, 1:w:n-w+1, 0:w-1, 0:w-1)
I = (X + dX) + (Y + dY - 1).*n
for k in 1:Base.div(n, w)
for l in 1:Base.div(n, w)
I[k, l, :, :] = transpose(I[k, l, :, :])
end
end
I
WARNING: Method definition ndgrid(AbstractArray{
64×64×4×4 Array{Int64,4}: [:, :, 1, 1] = 1 1025 2049 3073 4097 5121 … 60417 61441 62465 63489 64513 5 1029 2053 3077 4101 5125 60421 61445 62469 63493 64517 9 1033 2057 3081 4105 5129 60425 61449 62473 63497 64521 13 1037 2061 3085 4109 5133 60429 61453 62477 63501 64525 17 1041 2065 3089 4113 5137 60433 61457 62481 63505 64529 21 1045 2069 3093 4117 5141 … 60437 61461 62485 63509 64533 25 1049 2073 3097 4121 5145 60441 61465 62489 63513 64537 29 1053 2077 3101 4125 5149 60445 61469 62493 63517 64541 33 1057 2081 3105 4129 5153 60449 61473 62497 63521 64545 37 1061 2085 3109 4133 5157 60453 61477 62501 63525 64549 41 1065 2089 3113 4137 5161 … 60457 61481 62505 63529 64553 45 1069 2093 3117 4141 5165 60461 61485 62509 63533 64557 49 1073 2097 3121 4145 5169 60465 61489 62513 63537 64561 ⋮ ⋮ ⋱ ⋮ 209 1233 2257 3281 4305 5329 60625 61649 62673 63697 64721 213 1237 2261 3285 4309 5333 60629 61653 62677 63701 64725 217 1241 2265 3289 4313 5337 60633 61657 62681 63705 64729 221 1245 2269 3293 4317 5341 … 60637 61661 62685 63709 64733 225 1249 2273 3297 4321 5345 60641 61665 62689 63713 64737 229 1253 2277 3301 4325 5349 60645 61669 62693 63717 64741 233 1257 2281 3305 4329 5353 60649 61673 62697 63721 64745 237 1261 2285 3309 4333 5357 60653 61677 62701 63725 64749 241 1265 2289 3313 4337 5361 … 60657 61681 62705 63729 64753 245 1269 2293 3317 4341 5365 60661 61685 62709 63733 64757 249 1273 2297 3321 4345 5369 60665 61689 62713 63737 64761 253 1277 2301 3325 4349 5373 60669 61693 62717 63741 64765 [:, :, 2, 1] = 257 1281 2305 3329 4353 5377 … 60673 61697 62721 63745 64769 261 1285 2309 3333 4357 5381 60677 61701 62725 63749 64773 265 1289 2313 3337 4361 5385 60681 61705 62729 63753 64777 269 1293 2317 3341 4365 5389 60685 61709 62733 63757 64781 273 1297 2321 3345 4369 5393 60689 61713 62737 63761 64785 277 1301 2325 3349 4373 5397 … 60693 61717 62741 63765 64789 281 1305 2329 3353 4377 5401 60697 61721 62745 63769 64793 285 1309 2333 3357 4381 5405 60701 61725 62749 63773 64797 289 1313 2337 3361 4385 5409 60705 61729 62753 63777 64801 293 1317 2341 3365 4389 5413 60709 61733 62757 63781 64805 297 1321 2345 3369 4393 5417 … 60713 61737 62761 63785 64809 301 1325 2349 3373 4397 5421 60717 61741 62765 63789 64813 305 1329 2353 3377 4401 5425 60721 61745 62769 63793 64817 ⋮ ⋮ ⋱ ⋮ 465 1489 2513 3537 4561 5585 60881 61905 62929 63953 64977 469 1493 2517 3541 4565 5589 60885 61909 62933 63957 64981 473 1497 2521 3545 4569 5593 60889 61913 62937 63961 64985 477 1501 2525 3549 4573 5597 … 60893 61917 62941 63965 64989 481 1505 2529 3553 4577 5601 60897 61921 62945 63969 64993 485 1509 2533 3557 4581 5605 60901 61925 62949 63973 64997 489 1513 2537 3561 4585 5609 60905 61929 62953 63977 65001 493 1517 2541 3565 4589 5613 60909 61933 62957 63981 65005 497 1521 2545 3569 4593 5617 … 60913 61937 62961 63985 65009 501 1525 2549 3573 4597 5621 60917 61941 62965 63989 65013 505 1529 2553 3577 4601 5625 60921 61945 62969 63993 65017 509 1533 2557 3581 4605 5629 60925 61949 62973 63997 65021 [:, :, 3, 1] = 513 1537 2561 3585 4609 5633 … 60929 61953 62977 64001 65025 517 1541 2565 3589 4613 5637 60933 61957 62981 64005 65029 521 1545 2569 3593 4617 5641 60937 61961 62985 64009 65033 525 1549 2573 3597 4621 5645 60941 61965 62989 64013 65037 529 1553 2577 3601 4625 5649 60945 61969 62993 64017 65041 533 1557 2581 3605 4629 5653 … 60949 61973 62997 64021 65045 537 1561 2585 3609 4633 5657 60953 61977 63001 64025 65049 541 1565 2589 3613 4637 5661 60957 61981 63005 64029 65053 545 1569 2593 3617 4641 5665 60961 61985 63009 64033 65057 549 1573 2597 3621 4645 5669 60965 61989 63013 64037 65061 553 1577 2601 3625 4649 5673 … 60969 61993 63017 64041 65065 557 1581 2605 3629 4653 5677 60973 61997 63021 64045 65069 561 1585 2609 3633 4657 5681 60977 62001 63025 64049 65073 ⋮ ⋮ ⋱ ⋮ 721 1745 2769 3793 4817 5841 61137 62161 63185 64209 65233 725 1749 2773 3797 4821 5845 61141 62165 63189 64213 65237 729 1753 2777 3801 4825 5849 61145 62169 63193 64217 65241 733 1757 2781 3805 4829 5853 … 61149 62173 63197 64221 65245 737 1761 2785 3809 4833 5857 61153 62177 63201 64225 65249 741 1765 2789 3813 4837 5861 61157 62181 63205 64229 65253 745 1769 2793 3817 4841 5865 61161 62185 63209 64233 65257 749 1773 2797 3821 4845 5869 61165 62189 63213 64237 65261 753 1777 2801 3825 4849 5873 … 61169 62193 63217 64241 65265 757 1781 2805 3829 4853 5877 61173 62197 63221 64245 65269 761 1785 2809 3833 4857 5881 61177 62201 63225 64249 65273 765 1789 2813 3837 4861 5885 61181 62205 63229 64253 65277 [:, :, 4, 1] = 769 1793 2817 3841 4865 5889 … 61185 62209 63233 64257 65281 773 1797 2821 3845 4869 5893 61189 62213 63237 64261 65285 777 1801 2825 3849 4873 5897 61193 62217 63241 64265 65289 781 1805 2829 3853 4877 5901 61197 62221 63245 64269 65293 785 1809 2833 3857 4881 5905 61201 62225 63249 64273 65297 789 1813 2837 3861 4885 5909 … 61205 62229 63253 64277 65301 793 1817 2841 3865 4889 5913 61209 62233 63257 64281 65305 797 1821 2845 3869 4893 5917 61213 62237 63261 64285 65309 801 1825 2849 3873 4897 5921 61217 62241 63265 64289 65313 805 1829 2853 3877 4901 5925 61221 62245 63269 64293 65317 809 1833 2857 3881 4905 5929 … 61225 62249 63273 64297 65321 813 1837 2861 3885 4909 5933 61229 62253 63277 64301 65325 817 1841 2865 3889 4913 5937 61233 62257 63281 64305 65329 ⋮ ⋮ ⋱ ⋮ 977 2001 3025 4049 5073 6097 61393 62417 63441 64465 65489 981 2005 3029 4053 5077 6101 61397 62421 63445 64469 65493 985 2009 3033 4057 5081 6105 61401 62425 63449 64473 65497 989 2013 3037 4061 5085 6109 … 61405 62429 63453 64477 65501 993 2017 3041 4065 5089 6113 61409 62433 63457 64481 65505 997 2021 3045 4069 5093 6117 61413 62437 63461 64485 65509 1001 2025 3049 4073 5097 6121 61417 62441 63465 64489 65513 1005 2029 3053 4077 5101 6125 61421 62445 63469 64493 65517 1009 2033 3057 4081 5105 6129 … 61425 62449 63473 64497 65521 1013 2037 3061 4085 5109 6133 61429 62453 63477 64501 65525 1017 2041 3065 4089 5113 6137 61433 62457 63481 64505 65529 1021 2045 3069 4093 5117 6141 61437 62461 63485 64509 65533 [:, :, 1, 2] = 2 1026 2050 3074 4098 5122 … 60418 61442 62466 63490 64514 6 1030 2054 3078 4102 5126 60422 61446 62470 63494 64518 10 1034 2058 3082 4106 5130 60426 61450 62474 63498 64522 14 1038 2062 3086 4110 5134 60430 61454 62478 63502 64526 18 1042 2066 3090 4114 5138 60434 61458 62482 63506 64530 22 1046 2070 3094 4118 5142 … 60438 61462 62486 63510 64534 26 1050 2074 3098 4122 5146 60442 61466 62490 63514 64538 30 1054 2078 3102 4126 5150 60446 61470 62494 63518 64542 34 1058 2082 3106 4130 5154 60450 61474 62498 63522 64546 38 1062 2086 3110 4134 5158 60454 61478 62502 63526 64550 42 1066 2090 3114 4138 5162 … 60458 61482 62506 63530 64554 46 1070 2094 3118 4142 5166 60462 61486 62510 63534 64558 50 1074 2098 3122 4146 5170 60466 61490 62514 63538 64562 ⋮ ⋮ ⋱ ⋮ 210 1234 2258 3282 4306 5330 60626 61650 62674 63698 64722 214 1238 2262 3286 4310 5334 60630 61654 62678 63702 64726 218 1242 2266 3290 4314 5338 60634 61658 62682 63706 64730 222 1246 2270 3294 4318 5342 … 60638 61662 62686 63710 64734 226 1250 2274 3298 4322 5346 60642 61666 62690 63714 64738 230 1254 2278 3302 4326 5350 60646 61670 62694 63718 64742 234 1258 2282 3306 4330 5354 60650 61674 62698 63722 64746 238 1262 2286 3310 4334 5358 60654 61678 62702 63726 64750 242 1266 2290 3314 4338 5362 … 60658 61682 62706 63730 64754 246 1270 2294 3318 4342 5366 60662 61686 62710 63734 64758 250 1274 2298 3322 4346 5370 60666 61690 62714 63738 64762 254 1278 2302 3326 4350 5374 60670 61694 62718 63742 64766 [:, :, 2, 2] = 258 1282 2306 3330 4354 5378 … 60674 61698 62722 63746 64770 262 1286 2310 3334 4358 5382 60678 61702 62726 63750 64774 266 1290 2314 3338 4362 5386 60682 61706 62730 63754 64778 270 1294 2318 3342 4366 5390 60686 61710 62734 63758 64782 274 1298 2322 3346 4370 5394 60690 61714 62738 63762 64786 278 1302 2326 3350 4374 5398 … 60694 61718 62742 63766 64790 282 1306 2330 3354 4378 5402 60698 61722 62746 63770 64794 286 1310 2334 3358 4382 5406 60702 61726 62750 63774 64798 290 1314 2338 3362 4386 5410 60706 61730 62754 63778 64802 294 1318 2342 3366 4390 5414 60710 61734 62758 63782 64806 298 1322 2346 3370 4394 5418 … 60714 61738 62762 63786 64810 302 1326 2350 3374 4398 5422 60718 61742 62766 63790 64814 306 1330 2354 3378 4402 5426 60722 61746 62770 63794 64818 ⋮ ⋮ ⋱ ⋮ 466 1490 2514 3538 4562 5586 60882 61906 62930 63954 64978 470 1494 2518 3542 4566 5590 60886 61910 62934 63958 64982 474 1498 2522 3546 4570 5594 60890 61914 62938 63962 64986 478 1502 2526 3550 4574 5598 … 60894 61918 62942 63966 64990 482 1506 2530 3554 4578 5602 60898 61922 62946 63970 64994 486 1510 2534 3558 4582 5606 60902 61926 62950 63974 64998 490 1514 2538 3562 4586 5610 60906 61930 62954 63978 65002 494 1518 2542 3566 4590 5614 60910 61934 62958 63982 65006 498 1522 2546 3570 4594 5618 … 60914 61938 62962 63986 65010 502 1526 2550 3574 4598 5622 60918 61942 62966 63990 65014 506 1530 2554 3578 4602 5626 60922 61946 62970 63994 65018 510 1534 2558 3582 4606 5630 60926 61950 62974 63998 65022 [:, :, 3, 2] = 514 1538 2562 3586 4610 5634 … 60930 61954 62978 64002 65026 518 1542 2566 3590 4614 5638 60934 61958 62982 64006 65030 522 1546 2570 3594 4618 5642 60938 61962 62986 64010 65034 526 1550 2574 3598 4622 5646 60942 61966 62990 64014 65038 530 1554 2578 3602 4626 5650 60946 61970 62994 64018 65042 534 1558 2582 3606 4630 5654 … 60950 61974 62998 64022 65046 538 1562 2586 3610 4634 5658 60954 61978 63002 64026 65050 542 1566 2590 3614 4638 5662 60958 61982 63006 64030 65054 546 1570 2594 3618 4642 5666 60962 61986 63010 64034 65058 550 1574 2598 3622 4646 5670 60966 61990 63014 64038 65062 554 1578 2602 3626 4650 5674 … 60970 61994 63018 64042 65066 558 1582 2606 3630 4654 5678 60974 61998 63022 64046 65070 562 1586 2610 3634 4658 5682 60978 62002 63026 64050 65074 ⋮ ⋮ ⋱ ⋮ 722 1746 2770 3794 4818 5842 61138 62162 63186 64210 65234 726 1750 2774 3798 4822 5846 61142 62166 63190 64214 65238 730 1754 2778 3802 4826 5850 61146 62170 63194 64218 65242 734 1758 2782 3806 4830 5854 … 61150 62174 63198 64222 65246 738 1762 2786 3810 4834 5858 61154 62178 63202 64226 65250 742 1766 2790 3814 4838 5862 61158 62182 63206 64230 65254 746 1770 2794 3818 4842 5866 61162 62186 63210 64234 65258 750 1774 2798 3822 4846 5870 61166 62190 63214 64238 65262 754 1778 2802 3826 4850 5874 … 61170 62194 63218 64242 65266 758 1782 2806 3830 4854 5878 61174 62198 63222 64246 65270 762 1786 2810 3834 4858 5882 61178 62202 63226 64250 65274 766 1790 2814 3838 4862 5886 61182 62206 63230 64254 65278 [:, :, 4, 2] = 770 1794 2818 3842 4866 5890 … 61186 62210 63234 64258 65282 774 1798 2822 3846 4870 5894 61190 62214 63238 64262 65286 778 1802 2826 3850 4874 5898 61194 62218 63242 64266 65290 782 1806 2830 3854 4878 5902 61198 62222 63246 64270 65294 786 1810 2834 3858 4882 5906 61202 62226 63250 64274 65298 790 1814 2838 3862 4886 5910 … 61206 62230 63254 64278 65302 794 1818 2842 3866 4890 5914 61210 62234 63258 64282 65306 798 1822 2846 3870 4894 5918 61214 62238 63262 64286 65310 802 1826 2850 3874 4898 5922 61218 62242 63266 64290 65314 806 1830 2854 3878 4902 5926 61222 62246 63270 64294 65318 810 1834 2858 3882 4906 5930 … 61226 62250 63274 64298 65322 814 1838 2862 3886 4910 5934 61230 62254 63278 64302 65326 818 1842 2866 3890 4914 5938 61234 62258 63282 64306 65330 ⋮ ⋮ ⋱ ⋮ 978 2002 3026 4050 5074 6098 61394 62418 63442 64466 65490 982 2006 3030 4054 5078 6102 61398 62422 63446 64470 65494 986 2010 3034 4058 5082 6106 61402 62426 63450 64474 65498 990 2014 3038 4062 5086 6110 … 61406 62430 63454 64478 65502 994 2018 3042 4066 5090 6114 61410 62434 63458 64482 65506 998 2022 3046 4070 5094 6118 61414 62438 63462 64486 65510 1002 2026 3050 4074 5098 6122 61418 62442 63466 64490 65514 1006 2030 3054 4078 5102 6126 61422 62446 63470 64494 65518 1010 2034 3058 4082 5106 6130 … 61426 62450 63474 64498 65522 1014 2038 3062 4086 5110 6134 61430 62454 63478 64502 65526 1018 2042 3066 4090 5114 6138 61434 62458 63482 64506 65530 1022 2046 3070 4094 5118 6142 61438 62462 63486 64510 65534 [:, :, 1, 3] = 3 1027 2051 3075 4099 5123 … 60419 61443 62467 63491 64515 7 1031 2055 3079 4103 5127 60423 61447 62471 63495 64519 11 1035 2059 3083 4107 5131 60427 61451 62475 63499 64523 15 1039 2063 3087 4111 5135 60431 61455 62479 63503 64527 19 1043 2067 3091 4115 5139 60435 61459 62483 63507 64531 23 1047 2071 3095 4119 5143 … 60439 61463 62487 63511 64535 27 1051 2075 3099 4123 5147 60443 61467 62491 63515 64539 31 1055 2079 3103 4127 5151 60447 61471 62495 63519 64543 35 1059 2083 3107 4131 5155 60451 61475 62499 63523 64547 39 1063 2087 3111 4135 5159 60455 61479 62503 63527 64551 43 1067 2091 3115 4139 5163 … 60459 61483 62507 63531 64555 47 1071 2095 3119 4143 5167 60463 61487 62511 63535 64559 51 1075 2099 3123 4147 5171 60467 61491 62515 63539 64563 ⋮ ⋮ ⋱ ⋮ 211 1235 2259 3283 4307 5331 60627 61651 62675 63699 64723 215 1239 2263 3287 4311 5335 60631 61655 62679 63703 64727 219 1243 2267 3291 4315 5339 60635 61659 62683 63707 64731 223 1247 2271 3295 4319 5343 … 60639 61663 62687 63711 64735 227 1251 2275 3299 4323 5347 60643 61667 62691 63715 64739 231 1255 2279 3303 4327 5351 60647 61671 62695 63719 64743 235 1259 2283 3307 4331 5355 60651 61675 62699 63723 64747 239 1263 2287 3311 4335 5359 60655 61679 62703 63727 64751 243 1267 2291 3315 4339 5363 … 60659 61683 62707 63731 64755 247 1271 2295 3319 4343 5367 60663 61687 62711 63735 64759 251 1275 2299 3323 4347 5371 60667 61691 62715 63739 64763 255 1279 2303 3327 4351 5375 60671 61695 62719 63743 64767 [:, :, 2, 3] = 259 1283 2307 3331 4355 5379 … 60675 61699 62723 63747 64771 263 1287 2311 3335 4359 5383 60679 61703 62727 63751 64775 267 1291 2315 3339 4363 5387 60683 61707 62731 63755 64779 271 1295 2319 3343 4367 5391 60687 61711 62735 63759 64783 275 1299 2323 3347 4371 5395 60691 61715 62739 63763 64787 279 1303 2327 3351 4375 5399 … 60695 61719 62743 63767 64791 283 1307 2331 3355 4379 5403 60699 61723 62747 63771 64795 287 1311 2335 3359 4383 5407 60703 61727 62751 63775 64799 291 1315 2339 3363 4387 5411 60707 61731 62755 63779 64803 295 1319 2343 3367 4391 5415 60711 61735 62759 63783 64807 299 1323 2347 3371 4395 5419 … 60715 61739 62763 63787 64811 303 1327 2351 3375 4399 5423 60719 61743 62767 63791 64815 307 1331 2355 3379 4403 5427 60723 61747 62771 63795 64819 ⋮ ⋮ ⋱ ⋮ 467 1491 2515 3539 4563 5587 60883 61907 62931 63955 64979 471 1495 2519 3543 4567 5591 60887 61911 62935 63959 64983 475 1499 2523 3547 4571 5595 60891 61915 62939 63963 64987 479 1503 2527 3551 4575 5599 … 60895 61919 62943 63967 64991 483 1507 2531 3555 4579 5603 60899 61923 62947 63971 64995 487 1511 2535 3559 4583 5607 60903 61927 62951 63975 64999 491 1515 2539 3563 4587 5611 60907 61931 62955 63979 65003 495 1519 2543 3567 4591 5615 60911 61935 62959 63983 65007 499 1523 2547 3571 4595 5619 … 60915 61939 62963 63987 65011 503 1527 2551 3575 4599 5623 60919 61943 62967 63991 65015 507 1531 2555 3579 4603 5627 60923 61947 62971 63995 65019 511 1535 2559 3583 4607 5631 60927 61951 62975 63999 65023 [:, :, 3, 3] = 515 1539 2563 3587 4611 5635 … 60931 61955 62979 64003 65027 519 1543 2567 3591 4615 5639 60935 61959 62983 64007 65031 523 1547 2571 3595 4619 5643 60939 61963 62987 64011 65035 527 1551 2575 3599 4623 5647 60943 61967 62991 64015 65039 531 1555 2579 3603 4627 5651 60947 61971 62995 64019 65043 535 1559 2583 3607 4631 5655 … 60951 61975 62999 64023 65047 539 1563 2587 3611 4635 5659 60955 61979 63003 64027 65051 543 1567 2591 3615 4639 5663 60959 61983 63007 64031 65055 547 1571 2595 3619 4643 5667 60963 61987 63011 64035 65059 551 1575 2599 3623 4647 5671 60967 61991 63015 64039 65063 555 1579 2603 3627 4651 5675 … 60971 61995 63019 64043 65067 559 1583 2607 3631 4655 5679 60975 61999 63023 64047 65071 563 1587 2611 3635 4659 5683 60979 62003 63027 64051 65075 ⋮ ⋮ ⋱ ⋮ 723 1747 2771 3795 4819 5843 61139 62163 63187 64211 65235 727 1751 2775 3799 4823 5847 61143 62167 63191 64215 65239 731 1755 2779 3803 4827 5851 61147 62171 63195 64219 65243 735 1759 2783 3807 4831 5855 … 61151 62175 63199 64223 65247 739 1763 2787 3811 4835 5859 61155 62179 63203 64227 65251 743 1767 2791 3815 4839 5863 61159 62183 63207 64231 65255 747 1771 2795 3819 4843 5867 61163 62187 63211 64235 65259 751 1775 2799 3823 4847 5871 61167 62191 63215 64239 65263 755 1779 2803 3827 4851 5875 … 61171 62195 63219 64243 65267 759 1783 2807 3831 4855 5879 61175 62199 63223 64247 65271 763 1787 2811 3835 4859 5883 61179 62203 63227 64251 65275 767 1791 2815 3839 4863 5887 61183 62207 63231 64255 65279 [:, :, 4, 3] = 771 1795 2819 3843 4867 5891 … 61187 62211 63235 64259 65283 775 1799 2823 3847 4871 5895 61191 62215 63239 64263 65287 779 1803 2827 3851 4875 5899 61195 62219 63243 64267 65291 783 1807 2831 3855 4879 5903 61199 62223 63247 64271 65295 787 1811 2835 3859 4883 5907 61203 62227 63251 64275 65299 791 1815 2839 3863 4887 5911 … 61207 62231 63255 64279 65303 795 1819 2843 3867 4891 5915 61211 62235 63259 64283 65307 799 1823 2847 3871 4895 5919 61215 62239 63263 64287 65311 803 1827 2851 3875 4899 5923 61219 62243 63267 64291 65315 807 1831 2855 3879 4903 5927 61223 62247 63271 64295 65319 811 1835 2859 3883 4907 5931 … 61227 62251 63275 64299 65323 815 1839 2863 3887 4911 5935 61231 62255 63279 64303 65327 819 1843 2867 3891 4915 5939 61235 62259 63283 64307 65331 ⋮ ⋮ ⋱ ⋮ 979 2003 3027 4051 5075 6099 61395 62419 63443 64467 65491 983 2007 3031 4055 5079 6103 61399 62423 63447 64471 65495 987 2011 3035 4059 5083 6107 61403 62427 63451 64475 65499 991 2015 3039 4063 5087 6111 … 61407 62431 63455 64479 65503 995 2019 3043 4067 5091 6115 61411 62435 63459 64483 65507 999 2023 3047 4071 5095 6119 61415 62439 63463 64487 65511 1003 2027 3051 4075 5099 6123 61419 62443 63467 64491 65515 1007 2031 3055 4079 5103 6127 61423 62447 63471 64495 65519 1011 2035 3059 4083 5107 6131 … 61427 62451 63475 64499 65523 1015 2039 3063 4087 5111 6135 61431 62455 63479 64503 65527 1019 2043 3067 4091 5115 6139 61435 62459 63483 64507 65531 1023 2047 3071 4095 5119 6143 61439 62463 63487 64511 65535 [:, :, 1, 4] = 4 1028 2052 3076 4100 5124 … 60420 61444 62468 63492 64516 8 1032 2056 3080 4104 5128 60424 61448 62472 63496 64520 12 1036 2060 3084 4108 5132 60428 61452 62476 63500 64524 16 1040 2064 3088 4112 5136 60432 61456 62480 63504 64528 20 1044 2068 3092 4116 5140 60436 61460 62484 63508 64532 24 1048 2072 3096 4120 5144 … 60440 61464 62488 63512 64536 28 1052 2076 3100 4124 5148 60444 61468 62492 63516 64540 32 1056 2080 3104 4128 5152 60448 61472 62496 63520 64544 36 1060 2084 3108 4132 5156 60452 61476 62500 63524 64548 40 1064 2088 3112 4136 5160 60456 61480 62504 63528 64552 44 1068 2092 3116 4140 5164 … 60460 61484 62508 63532 64556 48 1072 2096 3120 4144 5168 60464 61488 62512 63536 64560 52 1076 2100 3124 4148 5172 60468 61492 62516 63540 64564 ⋮ ⋮ ⋱ ⋮ 212 1236 2260 3284 4308 5332 60628 61652 62676 63700 64724 216 1240 2264 3288 4312 5336 60632 61656 62680 63704 64728 220 1244 2268 3292 4316 5340 60636 61660 62684 63708 64732 224 1248 2272 3296 4320 5344 … 60640 61664 62688 63712 64736 228 1252 2276 3300 4324 5348 60644 61668 62692 63716 64740 232 1256 2280 3304 4328 5352 60648 61672 62696 63720 64744 236 1260 2284 3308 4332 5356 60652 61676 62700 63724 64748 240 1264 2288 3312 4336 5360 60656 61680 62704 63728 64752 244 1268 2292 3316 4340 5364 … 60660 61684 62708 63732 64756 248 1272 2296 3320 4344 5368 60664 61688 62712 63736 64760 252 1276 2300 3324 4348 5372 60668 61692 62716 63740 64764 256 1280 2304 3328 4352 5376 60672 61696 62720 63744 64768 [:, :, 2, 4] = 260 1284 2308 3332 4356 5380 … 60676 61700 62724 63748 64772 264 1288 2312 3336 4360 5384 60680 61704 62728 63752 64776 268 1292 2316 3340 4364 5388 60684 61708 62732 63756 64780 272 1296 2320 3344 4368 5392 60688 61712 62736 63760 64784 276 1300 2324 3348 4372 5396 60692 61716 62740 63764 64788 280 1304 2328 3352 4376 5400 … 60696 61720 62744 63768 64792 284 1308 2332 3356 4380 5404 60700 61724 62748 63772 64796 288 1312 2336 3360 4384 5408 60704 61728 62752 63776 64800 292 1316 2340 3364 4388 5412 60708 61732 62756 63780 64804 296 1320 2344 3368 4392 5416 60712 61736 62760 63784 64808 300 1324 2348 3372 4396 5420 … 60716 61740 62764 63788 64812 304 1328 2352 3376 4400 5424 60720 61744 62768 63792 64816 308 1332 2356 3380 4404 5428 60724 61748 62772 63796 64820 ⋮ ⋮ ⋱ ⋮ 468 1492 2516 3540 4564 5588 60884 61908 62932 63956 64980 472 1496 2520 3544 4568 5592 60888 61912 62936 63960 64984 476 1500 2524 3548 4572 5596 60892 61916 62940 63964 64988 480 1504 2528 3552 4576 5600 … 60896 61920 62944 63968 64992 484 1508 2532 3556 4580 5604 60900 61924 62948 63972 64996 488 1512 2536 3560 4584 5608 60904 61928 62952 63976 65000 492 1516 2540 3564 4588 5612 60908 61932 62956 63980 65004 496 1520 2544 3568 4592 5616 60912 61936 62960 63984 65008 500 1524 2548 3572 4596 5620 … 60916 61940 62964 63988 65012 504 1528 2552 3576 4600 5624 60920 61944 62968 63992 65016 508 1532 2556 3580 4604 5628 60924 61948 62972 63996 65020 512 1536 2560 3584 4608 5632 60928 61952 62976 64000 65024 [:, :, 3, 4] = 516 1540 2564 3588 4612 5636 … 60932 61956 62980 64004 65028 520 1544 2568 3592 4616 5640 60936 61960 62984 64008 65032 524 1548 2572 3596 4620 5644 60940 61964 62988 64012 65036 528 1552 2576 3600 4624 5648 60944 61968 62992 64016 65040 532 1556 2580 3604 4628 5652 60948 61972 62996 64020 65044 536 1560 2584 3608 4632 5656 … 60952 61976 63000 64024 65048 540 1564 2588 3612 4636 5660 60956 61980 63004 64028 65052 544 1568 2592 3616 4640 5664 60960 61984 63008 64032 65056 548 1572 2596 3620 4644 5668 60964 61988 63012 64036 65060 552 1576 2600 3624 4648 5672 60968 61992 63016 64040 65064 556 1580 2604 3628 4652 5676 … 60972 61996 63020 64044 65068 560 1584 2608 3632 4656 5680 60976 62000 63024 64048 65072 564 1588 2612 3636 4660 5684 60980 62004 63028 64052 65076 ⋮ ⋮ ⋱ ⋮ 724 1748 2772 3796 4820 5844 61140 62164 63188 64212 65236 728 1752 2776 3800 4824 5848 61144 62168 63192 64216 65240 732 1756 2780 3804 4828 5852 61148 62172 63196 64220 65244 736 1760 2784 3808 4832 5856 … 61152 62176 63200 64224 65248 740 1764 2788 3812 4836 5860 61156 62180 63204 64228 65252 744 1768 2792 3816 4840 5864 61160 62184 63208 64232 65256 748 1772 2796 3820 4844 5868 61164 62188 63212 64236 65260 752 1776 2800 3824 4848 5872 61168 62192 63216 64240 65264 756 1780 2804 3828 4852 5876 … 61172 62196 63220 64244 65268 760 1784 2808 3832 4856 5880 61176 62200 63224 64248 65272 764 1788 2812 3836 4860 5884 61180 62204 63228 64252 65276 768 1792 2816 3840 4864 5888 61184 62208 63232 64256 65280 [:, :, 4, 4] = 772 1796 2820 3844 4868 5892 … 61188 62212 63236 64260 65284 776 1800 2824 3848 4872 5896 61192 62216 63240 64264 65288 780 1804 2828 3852 4876 5900 61196 62220 63244 64268 65292 784 1808 2832 3856 4880 5904 61200 62224 63248 64272 65296 788 1812 2836 3860 4884 5908 61204 62228 63252 64276 65300 792 1816 2840 3864 4888 5912 … 61208 62232 63256 64280 65304 796 1820 2844 3868 4892 5916 61212 62236 63260 64284 65308 800 1824 2848 3872 4896 5920 61216 62240 63264 64288 65312 804 1828 2852 3876 4900 5924 61220 62244 63268 64292 65316 808 1832 2856 3880 4904 5928 61224 62248 63272 64296 65320 812 1836 2860 3884 4908 5932 … 61228 62252 63276 64300 65324 816 1840 2864 3888 4912 5936 61232 62256 63280 64304 65328 820 1844 2868 3892 4916 5940 61236 62260 63284 64308 65332 ⋮ ⋮ ⋱ ⋮ 980 2004 3028 4052 5076 6100 61396 62420 63444 64468 65492 984 2008 3032 4056 5080 6104 61400 62424 63448 64472 65496 988 2012 3036 4060 5084 6108 61404 62428 63452 64476 65500 992 2016 3040 4064 5088 6112 … 61408 62432 63456 64480 65504 996 2020 3044 4068 5092 6116 61412 62436 63460 64484 65508 1000 2024 3048 4072 5096 6120 61416 62440 63464 64488 65512 1004 2028 3052 4076 5100 6124 61420 62444 63468 64492 65516 1008 2032 3056 4080 5104 6128 61424 62448 63472 64496 65520 1012 2036 3060 4084 5108 6132 … 61428 62452 63476 64500 65524 1016 2040 3064 4088 5112 6136 61432 62456 63480 64504 65528 1020 2044 3068 4092 5116 6140 61436 62460 63484 64508 65532 1024 2048 3072 4096 5120 6144 61440 62464 63488 64512 65536
T<:Any, 1}) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:3 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:3. WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:6 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:6. WARNING: Method definition ndgrid_fill(Any, Any, Any, Any) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:13 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:13. WARNING: Method definition ndgrid(AbstractArray{#T<:Any, 1}...) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:19 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:19. WARNING: Method definition meshgrid(AbstractArray{T<:Any, 1}) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:33 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:33. WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:36 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:36. WARNING: Method definition meshgrid(AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}, AbstractArray{#T<:Any, 1}) in module Main at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:44 overwritten at C:\Users\Ayman\.julia\v0.5\ndgrid.jl:44.
Block extraction operator. It returns the set $ \{x_b\}_{b \in B} $ of block-partitioned coefficients.
block = x -> x[I]
(::#93) (generic function with 1 method)
Block reconstruction operator.
function assign(M,I,H)
M_temp = M
M_temp[I] = H
return reshape(M_temp, n, n)
end
iblock = H -> assign(zeros(n, n), I, H)
(::#95) (generic function with 1 method)
WARNING: Method definition assign(Any, Any, Any) in module Main at In[100]:2 overwritten at In[121]:2.
Check that block extraction / reconstruction gives perfect reconstruction.
print(string("Should be 0:", string(norm(f - iblock(block(f))))))
Should be 0:0.0
Compute the average energy of each block, and duplicate.
function energy(H)
H_tmp = copy(H)
for i in 1:Base.div(n, w)
for j in 1:Base.div(n, w)
H_tmp[i, j, :, :] = sqrt(mean(H_tmp[i, j, :, :].^2))#*np.ones([1,1])
end
end
return H_tmp
end
energy (generic function with 1 method)
Block thresholding operator.
Thresh = (H, T) -> psi(energy(H), T).*H
ThreshBlock = (x, T) -> iblock(Thresh(block(x), T))
(::#99) (generic function with 1 method)
Exercise 2
Test the effect of block thresholding on the image $f_0$ itself, for increasing value of $T$. Of course directly thresholding the image has no interest, this is just to vizualize the effect.
# run -i nt_solutions/denoisingwav_4_block/exo2
include("Exos\\denoisingwav_4_block\\exo2.jl")
## Insert your code here.
Wavelet coefficients of natural images are not independant one from each other. One can thus improve the denoising results by thresholding block of coefficients togethers. Block thresholding is only efficient when used as a soft thresholder. Here we use a Stein soft thresholder.
Display the thresholded coefficients for a threshold value $T$ proportional to the noise level $\si$.
T = 1.25*sigma
figure(figsize = (10, 10))
plot_wavelet(ThreshBlock(wav(f), T), Jmin)
show()
Define the wavelet block thresholding operator.
ThreshWav = (f, T) -> iwav(ThreshBlock(wav(f), T))
(::#103) (generic function with 1 method)
Test the thresholding.
figure(figsize = (5, 5))
imageplot(clamP(ThreshWav(f, T)))
BoundsError: attempt to access 32×32 Array{Float64,2} at index [1:512,Colon()] in throw_boundserror(::Array{Float64,2}, ::Tuple{UnitRange{Int64},Colon}) 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 lifting_step(::Array{Float64,2}, ::Array{Float64,1}, ::Int64) at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:179 in perform_wavelet_transf(::Array{Float64,2}, ::Int64, ::Int64, ::String, ::Int64, ::Int64) at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:90 in (::##3#4)(::Array{Float64,2}) at .\In[10]:3 in (::##103#104)(::Array{Float64,2}, ::Float64) at .\In[140]:1
Exercise 3
Display the evolution of the denoising SNR when $T$ varies. Store the optimal denoising result in $f_{Block}$.
run -i nt_solutions/denoisingwav_4_block/exo3
## Insert your code here.
Display the result.
plt.figure(figsize=(5,5))
imageplot(clamp(fBlock), "SNR = %.1f dB" %snr(f0, fBlock))
Block thresholding can also be applied to a translation invariant wavelet transform. It gives state of the art denoising results.
Shortcuts for the foward and backward translation invariant wavelet transforms.
wav = f -> perform_wavelet_transf(f, Jmin, + 1, "9-7", 0, 1)
iwav = fw -> perform_wavelet_transf(fw, Jmin, -1,"9-7", 0, 1)
(::#115) (generic function with 1 method)
Foward wavelet transform.
fw = wav(f)
n = 256
#np.shape(fw)[0]
BoundsError: attempt to access 256×13×1 Array{Float64,3} at index [Colon(),[2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,255],Colon()] in throw_boundserror(::Array{Float64,3}, ::Tuple{Colon,Array{Int64,1},Colon}) at .\abstractarray.jl:355 in checkbounds at .\abstractarray.jl:284 [inlined] in _getindex at .\multidimensional.jl:270 [inlined] in getindex(::Array{Float64,3}, ::Colon, ::Array{Int64,1}, ::Colon) at .\abstractarray.jl:752 in lifting_step_ti(::Array{Float64,2}, ::Array{Float64,1}, ::Int64, ::Int64) at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:227 in perform_wavelet_transf(::Array{Float64,2}, ::Int64, ::Int64, ::String, ::Int64, ::Int64) at C:\Users\Ayman\.julia\v0.5\NtToolBox\src\perform_wavelet_transf.jl:118 in (::##113#114)(::Array{Float64,2}) at .\In[146]:1
Compute indexing of the blocks.
[X,J,Y,dX,dY] = np.meshgrid(np.arange(1,n-w+2,w),np.arange(1,np.shape(fw)[0]+1),np.arange(1,n-w+2,w),np.arange(0,w),np.arange(0,w))
I = (X + dX-1) + (Y + dY-1)*n + (J-1)*n**2
for k in range(n//w):
for l in range(n//w):
for m in range(np.shape(fw)[0]):
I[m][k][l] = np.transpose(I[m][k][l])
Forward and backward extraction operators.
block = lambda x : np.ravel(x)[I]
def assign(M,I,H):
M_temp = M
np.ravel(M_temp)[I] = H
return np.reshape(M_temp,(np.shape(fw)[0],n,n))
iblock = lambda H : assign(np.zeros([np.shape(fw)[0],n,n]), I, H)
Compute the average energy of each block, and duplicate.
def energy(H):
H_tmp = np.copy(H)
for i in range(n//w):
for j in range(n//w):
for k in range(np.shape(fw)[0]):
H_tmp[k][i][j] = np.sqrt(np.mean(H_tmp[k][i][j]**2))
return H_tmp
Block thresholding operator.
Thresh = lambda H,T : psi(energy(H), T)*H
ThreshBlock = lambda x,T : iblock(Thresh(block(x), T))
Define the wavelet block thresholding operator.
ThreshWav = lambda f,T : iwav(ThreshBlock(wav(f), T))
Test the thresholding.
T = 1.25*sigma
plt.figure(figsize = (5,5))
imageplot(clamp(ThreshWav(f, T)))
Exercise 4
Display the evolution of the denoising SNR when $T$ varies. Store the optimal denoising result in $f_{TI}$.
run -i nt_solutions/denoisingwav_4_block/exo4
## Insert your code here.
Display the result.
plt.figure(figsize = (5,5))
imageplot(clamp(fTI), "SNR = %.1f dB" %snr(f0, fTI))