-
Notifications
You must be signed in to change notification settings - Fork 0
/
HW03.jl
2251 lines (1808 loc) · 71.5 KB
/
HW03.jl
1
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
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
### A Pluto.jl notebook ###
# v0.19.14
using Markdown
using InteractiveUtils
# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
quote
local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
local el = $(esc(element))
global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
el
end
end
# ╔═╡ 77e6389e-3491-45b6-9f76-459b15a8922d
using TestImages, ImageShow, Random, PlutoTest, Plots, FourierTools, Statistics, Colors
# ╔═╡ ee4b6b2d-119e-4ac3-a1f5-500ff045bfc5
using FFTW, PlutoUI
# ╔═╡ c7150c2c-afeb-4158-b381-3341bb678fd5
using Noise
# ╔═╡ 31c31237-c7bd-4fdb-8f7f-219d323f80e7
"""
my_show(arr::AbstractArray{<:Real}; set_one=false, set_zero=false)
Displays a real valued array . Brightness encodes magnitude.
Works within Jupyter and Pluto.
## Keyword args
* `set_one=false` divides by the maximum to set maximum to 1
* `set_zero=false` subtracts the minimum to set minimum to 1
"""
function gray_show(arr::AbstractArray{<:Real}; set_one=true, set_zero=false)
arr = set_zero ? arr .- minimum(arr) : arr
arr = set_one ? arr ./ maximum(arr) : arr
Gray.(arr)
end
# ╔═╡ 7148e705-8b90-46ea-8763-778897daf9a7
md"# Homework 03
In this homework we cover several topics like image correction, Fourier transforms, Fourier Shifting and Sampling.
"
# ╔═╡ 83d9c252-2a8d-47a3-b247-a62b6a175201
md"## 1. Flat Fielding
As discussed in the lecture, raw images might be affected by dirty optics or sensor errors.
In this exercise we want to apply such a correction to simulated images.
"
# ╔═╡ 39c623c7-6063-4360-a2b7-95ea42c5262d
begin
img = Float32.(testimage("mandril_gray.tif"))
gray_show(img)
end
# ╔═╡ eacd8fdf-971c-46ca-914f-f73f65e28ca8
md"### Simulate Deteriorated Image
In the first step, we degrade our image with some artifacts.
We add an offset to our image which is caused by sensor artifacts.
Second, we add some dirty spots to the image which are for example caused by dust on the sensor.
"
# ╔═╡ 25dde5cb-9ed6-498e-b511-50ed0640d0c0
"""
simulate_bad_microscope(img)
The ideal `img` is degraded by artifacts.
"""
function simulate_bad_microscope(img::AbstractArray{T}) where T
# generate same random numbers!
rng = MersenneTwister(2021);
# some weird dark field slope offset
dark = zero.(img)
dark .+= range(0,0.3, length=size(img, 1))
xs = 1:size(img, 1)
ff = one.(img)
for i = 1:20 # emulates some dust spots as modulations to the flatfield ff
x, y = rand(rng, xs), rand(rng, xs)
sca = abs(randn(rng) .* T(0.2))
r = sqrt.((xs .- x).^2 .+ (xs' .- y).^2)
ff .-= T(0.2) .* sinc.(r .* sca)
end
return abs.(dark .+ ff .* img)
end
# ╔═╡ a6c024a4-7cb2-4f8d-be24-8727d75679da
img_dirty = simulate_bad_microscope(img);
# ╔═╡ d5c127ee-1ecf-4ed0-9cd9-6fcdefa95847
gray_show(img_dirty)
# ╔═╡ 6e95cbcb-0d89-4eca-8599-b266076708a7
md"## 1.1 Task
Extract both a flat field and a dark field image with `simulate_bad_microscope`.
Use the function `simulate_bad_microscope` to obtain a measured flat field and measured dark field.
"
# ╔═╡ 910cb7fa-b1fc-4106-acc4-3bf56edfb83d
flat_field = similar(img); # TODO
# ╔═╡ cbebb311-c521-4fa6-a909-1a4df8da4d2d
gray_show(flat_field)
# ╔═╡ 8bb737cb-2d59-49ba-9b89-afef6ca7bbf8
dark_field = similar(img); # TODO
# ╔═╡ 52bc8b3d-542e-4a35-b38c-35840f4075d4
gray_show(dark_field)
# ╔═╡ e3448ffe-5fbf-4bed-bf03-6eb9031e6c12
"""
correct_image(img_dirty, dark_field, flat_field)
Correct a dirty image `img_dirty` by providing the `dark_field` and `flat_field`.
Conceptually we divide by the `flat_field` since the measured `flat_field`
describes the deviations from an ideal white image.
The `dark_field` has to be subtracted before since that only contains sensor artifacts.
"""
function correct_image(img_dirty, dark_field, flat_field)
# TODO
return img_dirty
end
# ╔═╡ d9308666-efaa-4fcc-a083-f3d12dbab1e2
img_fixed = similar(img_dirty); # todo
# ╔═╡ 493b6bb4-1dfc-4a82-9360-40944b29826c
gray_show(img_fixed)
# ╔═╡ c80a2ce8-6620-4de1-bb78-5af97caf466d
md"### 1.1 Test"
# ╔═╡ e2e470e4-ff73-4dd5-973c-14e1535b01ea
PlutoTest.@test img_fixed ≈ img
# ╔═╡ 0a957f4a-a975-4087-b458-3236940eda43
PlutoTest.@test img ≈ correct_image(img, zero.(img), one.(img))
# ╔═╡ af9866f4-2b03-4dc4-9f8a-e02ae46e2f9c
md"# 2 Sampling
If a signal is not properly sampled, we can see a lot of artifacts which are called aliasing.
Try to determine experimentally (via the slider) which is the correct value for Nyquist sampling.
Afterwards, calculate the correct sampling in the variable $N$!
Don't state the pure number but rather a short equation which Julia calculates (something like `1 * 2 * 3`).
Do you see a difference between your expected and your calculated value?
"
# ╔═╡ 7a481af3-4ec8-4a58-bb84-e9d3d5381570
md"N=$(@bind N Slider(2:300, show_value=true))"
# ╔═╡ 4006ccf7-df8e-440c-8b3c-af525bb50ed4
begin
xs = range(0, 8, length=N+1)[begin:end-1]
f(x) = sin(x * 2π * 5)
end
# ╔═╡ c871cdbd-b560-46f2-918c-45d24d2f9b75
plot(xs, f.(xs))
# ╔═╡ f6b8cd92-c1ab-445b-8b4f-04f521b72cad
N_correct = 2 # TODO
# ╔═╡ 51b2d7be-fe94-44e9-8e67-73aafb636ef1
begin
xs_correct = range(0, 8, length=round(Int, N_correct + 1))[begin:end-1]
plot(xs_correct, f.(xs_correct))
end
# ╔═╡ a19c8850-9a7d-4e26-908b-a067a9006e84
md"No tests, since they would reveal the answer."
# ╔═╡ b0c67bd6-b5ef-4359-8727-1a3bfefcf759
md"# 3 FFTs
Another very important topics are FFTs!
The standard library in many languages is based on FFTW.
In Julia FFTW.jl is the recommended library.
The discrete Fourier transform is given by
$$X_k = \sum_{n=0}^{N-1} x_n \exp\left(-i2 \pi \frac{kn}{N} \right)$$
and the inverse by
$$x_n = \frac1{N}\sum_{k=0}^{N-1} X_k \exp\left(i2 \pi \frac{kn}{N} \right)$$.
The FFT algorithm evaluates those equations very efficiently ($\mathcal O(N \cdot \log N)$) by avoiding the explicit sum ($\mathcal O(N^2)$).
Further, the FFT calculates the frequencies in an (maybe) unexpected format
* Data: $[x_1, x_2, ..., x_N]$
* FFT of Data: $[f_0, f_1, ..., f_{N/2}, f_{-N/2}, f_{-N/2+1}, ...,f_{-1}]$
* (for real data negative and positive frequencies are the complex conjugate)
Sometimes we want to have the frequencies centered, which means a `fftshift` operations is needed
* FFT of Data followed by fftshift: $[f_{-N/2}, f_{-N/2+1},..., f_{-1}, f_0, f_1,..., f_{N/2-1}, f_{N/2}]$
"
# ╔═╡ 7276f9e0-80b8-4794-9125-fca43246f818
"""
my_fft(data::Vector{T})::Vector{Complex{T}}
Calculate the FFT of `data` according to the equation given above.
Note, the output should be a complex Vector!
"""
function my_fft(data::Vector{T})::Vector{Complex{T}} where T
out = zeros(Complex{T}, size(data))
# TODO
return out::Vector{Complex{T}}
end
# ╔═╡ e895cff6-4018-495f-9a0e-1e56ddcfd17d
arr = [1f0,2,3,4]
# ╔═╡ 14d68904-4de5-4f1d-81f6-dae8bf50b749
my_fft(arr)
# ╔═╡ 7307dd33-cb79-48c0-9e73-78ed70e9c628
fft(arr)
# ╔═╡ 393d03e2-461d-480a-9719-f33a4eef3a6f
md"### 3.1 Test"
# ╔═╡ df12fd86-f29d-4926-b360-daf4cb40af86
begin
arr_r = randn((13,))
PlutoTest.@test fft(arr_r) ≈ my_fft(arr_r)
end
# ╔═╡ 37eedc2e-27bd-4d1c-a83b-6ec14aba55ab
begin
arr_r2 = randn((14,))
PlutoTest.@test fft(arr_r2) ≈ my_fft(arr_r2)
end
# ╔═╡ 5313dbed-bfd7-4bb8-b39c-89c55e50d1bb
md"## 3.2 Task
Calculate the FFT of the following dataset, such that center frequency is in the center!
Use `fft` and `fftshift`.
"
# ╔═╡ e1eac3cd-e6ce-4743-9119-5ee0c69b8e00
begin
data_odd = [1.7142857142857142, 1.643544850921574, 1.0471849694228637, 1.383330948397344, 1.8300917955219322, 1.2951185346328968, 1.086443186817675]
data_even = [-2.5, 1.5, -2.5, 3.5]
end
# ╔═╡ 481961ac-ad7c-4768-a7aa-eed34143226b
ffts(x) = similar(x)# TODO
# ╔═╡ fc3bb52b-f9a2-42c0-b34b-347d0a137e2d
ffts(data_even)
# ╔═╡ a55851b6-3f8d-4af3-8f09-a3777c56df1b
ffts(data_odd)
# ╔═╡ e46d67be-d974-4f07-8b3e-02fff3e749d0
md"## 3.3 Task
Now do the inverse operation such that `data ≈ iffts(ffts(x))`
Check out `ifftshift`.
"
# ╔═╡ deb9eaf0-ee89-49c5-90b8-570d1a88f8fc
iffts(x) = similar(x) # TODO
# ╔═╡ 2875999a-a61d-491a-8a42-0552a87d3c6a
md"""
### 3.2 and 3.3 Test
"""
# ╔═╡ 333a253e-ce8b-477a-9db2-f99c034ae048
PlutoTest.@test real(ffts(iffts(data_even))) ≈ data_even
# ╔═╡ e81990c8-2194-472c-a102-b0d010f3a266
PlutoTest.@test real(ffts(iffts(data_odd))) ≈ data_odd
# ╔═╡ f3fc2ea5-033c-4789-8ce6-0c6b15b607b7
PlutoTest.@test real(iffts(ffts(data_odd))) ≈ data_odd
# ╔═╡ 90596b8a-53d4-457f-9e93-522e6c5dc35d
PlutoTest.@test real(iffts(ffts(data_even))) ≈ data_even
# ╔═╡ 4e54ddc6-b9cd-472f-b496-18055ee5b667
md"## 3.4 Task
Try to reproduce this $\cos$ curve.
The reason it looks so _sharp_ is that we only sample slightly above the Nyquist limit.
However, you are only allowed to use `fft`, `fftshift`, `ifftshift`, `ifft` and array manipulation (like `x[34] = 1723.12`).
And you can also take the `real` part at the end.
Please explain with some #comments what you did.
"
# ╔═╡ 9844ec03-b169-40a0-ad56-1c07b11bb886
begin
x = range(0, 4, length=21)[begin:end-1]
y = cos.(x .* 2π)
plot(x, y, mark="-*")
end
# ╔═╡ 49c39cea-895e-429a-909c-d1438d0dfd00
function reproduce_cos()
# still wrong
guess = zeros(20)
guess[3] = 20
return real(ift(guess))
end
# ╔═╡ 5c4f40df-cde3-4140-ace2-6c0e1b60ebcf
begin
plot(x, y, mark="-*")
plot!(x, reproduce_cos(), mark="-*")
end
# ╔═╡ 79578e4d-7527-459b-8915-ee5cb9ccad24
md"### 3.4 Test"
# ╔═╡ a74c4034-fa19-4360-9f5e-191d5ba4001a
PlutoTest.@test real(reproduce_cos()) ≈ y
# ╔═╡ 8ce249a3-0881-49e5-8ffd-5df21bda50f4
md"## 3.5 Fourier Shift Theorem
The discrete Fourier transform is given by
$$X_k = \sum_{n=0}^{N-1} x_n \exp\left(-i2 \pi \frac{kn}{N} \right)$$
and the inverse by
$$x_n = \frac1{N}\sum_{k=0}^{N-1} X_k \exp\left(i2 \pi \frac{kn}{N} \right)$$.
"
# ╔═╡ 13131381-11b3-4db4-9351-fb7264c94ea4
md"## 3.4 Task
Proof the Fourier shift theorem (with LaTeX) which states
$x_{n+\Delta n} = \mathcal{F}^{-1}\left[ X_{k} \cdot f(k, \Delta n)\right]$
Derive how the function $f(k, \Delta n)$ looks like!
hint: by clicking on the eye symbol next to this task, you can see how to embed LaTeX code in markdown comments.
"
# ╔═╡ bc3ae07f-3877-480f-bd0b-633ce21172e3
md"### Derivation
# TODO
"
# ╔═╡ 8bf6ba5b-1d06-41de-b22a-b58c1eaf2c83
"""
shift(x::AbstractArray{T, N}, Δn::Number)
Shifts a (real) array via the Fourier shift theorem by the amount of
`Δn`.
If `Δn ∈ ℕ` then `circshift(x, -Δn) ≈ shift(x, Δn).`
Otherwise a sub-pixel shift is obtained which is equivalent
to a sinc interpolation.
**Hint**: Use `fftfreq` to obtain the frequencies in Fourier space. Otherwise you might run into issues.
## Examples
```julia
julia> shift([1.0, 2.0, 3.0, 4.0], 1)
4-element Vector{Float64}:
2.0
3.0
4.0
1.0
julia> shift([1.0, 2.0, 3.0, 4.0], -1)
4-element Vector{Float64}:
4.0
1.0
2.0
3.0
julia> shift([1.0, 2.0, 3.0, 4.0], -0.5)
4-element Vector{Float64}:
2.5
1.085786437626905
2.5
3.914213562373095
```
"""
function shift(x::AbstractArray{T, N}, Δn::Number) where {T, N}
# TODO
return similar(x)
end
# ╔═╡ a21d8b13-20db-4df1-a2c7-451f8942466b
begin
x2 = 0:0.1:3
y2 = sin.(x2.*3) .+ cos.(x2 .* 4).^2
end;
# ╔═╡ 420d1b15-8dc9-4c9a-b91f-9cc4b8101b92
md"
$\Delta n$=$(@bind Δn Slider(0:0.01:100, show_value=true))"
# ╔═╡ 54ebf4c4-6692-4a40-a83e-c548440b3b85
begin
plot(shift(y2, Δn), label="FFT based shift", mark="-*")
plot!(circshift(y2, round(Int, .-Δn)), label="rounded to integer with circshift", mark="-*")
end
# ╔═╡ bbe52a98-9811-4c3a-9b4e-b90a1f5c163f
md"### 3.5 Test"
# ╔═╡ 73fd7360-ffdf-4870-a301-e39a982e0517
arr4 = randn((37,));
# ╔═╡ 37a815a0-c6cd-4ee8-9f46-ee4aa5259690
PlutoTest.@test circshift(arr4, -12) ≈ shift(arr4, 12)
# ╔═╡ 3c09615a-e9c4-462b-9f23-5d5676aa29b7
PlutoTest.@test FourierTools.shift(arr4, -13.1) ≈ shift(arr4, 13.1)
# ╔═╡ d3f2648a-a736-4812-b721-70298151a1a7
md"# 4 Noise Beyond Bandlimit
As seen in the lecture, noise can move information even beyond
the bandlimit of the optical system.
In this task, we want to demonstrate this behaviour with a simulation.
For that we need a few helper functions, try to fill the missing lines.
## 4.1 rr2 Function
In the first step, we create a rr2 function.
"
# ╔═╡ 9ef7592b-92f9-48bb-aa3e-9710bf251379
begin
"""
rr2(T, s, center=s .÷ 2 + 1)
Returns a 2D array which stores the squared distance to the center.
## Examples
```julia
julia> rr2(Float32, (5,5))
5×5 Matrix{Float32}:
8.0 5.0 4.0 5.0 8.0
5.0 2.0 1.0 2.0 5.0
4.0 1.0 0.0 1.0 4.0
5.0 2.0 1.0 2.0 5.0
8.0 5.0 4.0 5.0 8.0
julia> rr2((4,4))
4×4 Matrix{Float64}:
8.0 5.0 4.0 5.0
5.0 2.0 1.0 2.0
4.0 1.0 0.0 1.0
5.0 2.0 1.0 2.0
julia> rr2((3,3), (3,3))
3×3 Matrix{Float64}:
8.0 5.0 4.0
5.0 2.0 1.0
4.0 1.0 0.0
```
"""
function rr2(T::DataType, s, center=s .÷ 2 .+ 1)
# TODO
return zeros(s)
end
function rr2(s, center=s .÷ 2 .+ 1)
rr2(Float64, s, center)
end
end
# ╔═╡ b62519f5-26dd-49ca-afa4-2c72de7ef8a6
rr2((4,4))
# ╔═╡ bd45654d-15e0-4717-962b-d5a302de4a9b
PlutoTest.@test rr2(Float32, (4, 3), (2.5, 3)) == Float32[6.25 3.25 2.25; 4.25 1.25 0.25; 4.25 1.25 0.25; 6.25 3.25 2.25]
# ╔═╡ c6ed3889-0d44-4180-9f81-3095055c6e54
PlutoTest.@test rr2((3, 3), (3, 3)) == [8.0 5.0 4.0; 5.0 2.0 1.0; 4.0 1.0 0.0]
# ╔═╡ 5a16e94e-b01b-4897-9631-41d345229be2
md"## 4.2 circ Function
Having a function which creates the radial distance, we now create a function which returns an array which is 1 inside an circle of radius `r` and 0 outside of it.
"
# ╔═╡ 06a666f9-c99e-41e5-b703-a5e826721d04
begin
"""
circ(size, radius)
`size` the size of the resulting array
and `radius` the radius of the circle
"""
circ(size, radius) = ones(size) # TODO
end
# ╔═╡ 592e0ed4-bd97-4bb6-a826-b1e2f5ddbd59
circ((5,5), 2)
# ╔═╡ 205a744b-e8a1-4acc-ba6b-7a7307f9a028
md"## 4.2 Test"
# ╔═╡ daf4e474-0db4-4865-859d-96c31e7947f0
PlutoTest.@test circ((10, 2), 2.76) ≈ Bool[0 0; 0 0; 0 0; 1 1; 1 1; 1 1; 1 1; 1 1; 0 0; 0 0]
# ╔═╡ 7fb8fedf-5684-4874-8f17-d01debaff477
PlutoTest.@test circ((4, 5), 2.1) ≈ Bool[0 0 1 0 0; 0 1 1 1 0; 1 1 1 1 1; 0 1 1 1 0]
# ╔═╡ 9b9fb85b-13d0-43a2-ac5f-713593117b96
md" ## 4.3 Frequency Filter
Now we combine all methods together such that we frequency filter
our image.
Procedure:
* Transform image to Fourier space
* take only the inner part of a circle a `radius` of the Fourier space
* go to real space and take real part
"
# ╔═╡ 2a48072b-fc5a-4044-9aa9-5661483ffbfc
function frequency_filter(img, radius)
# TODO
return similar(img)
end
# ╔═╡ 643bf329-cce9-4c61-931d-65455b59357d
md"r=$(@bind r Slider(1:512, show_value=true))"
# ╔═╡ 2cafea2d-b4b6-4ede-b899-720a71034fda
gray_show(frequency_filter(img, r))
# ╔═╡ 7467608a-ce6d-4f29-8cd8-d766395f7129
md"## 4.4 Noise Breaks Resolution Limit?
We demonstrate now, that noise beyond the frequency limit still contains information
### Procedure
* First frequency filter our image such that it definetely bandlimited
* Apply Poissonian and Gaussian noise
* Now do the exact opposite of the above low-pass filter and take only the high frequency content!
"
# ╔═╡ 6d7fd369-1f6c-487d-b399-c8cbbfede916
function simulate(img, noise_function, radius)
# TODO
return similar(img)
end
# ╔═╡ e8bd9cb9-bd20-4619-93da-20e7829a61b4
# anonymous noise function we are going to pass to `simulate`
begin
p(x) = poisson(x, 100000)
g(x) = add_gauss(x, 0.1)
end
# ╔═╡ 8a479712-7290-4c00-af7c-44e11ec81d0d
md"r=$(@bind r2 Slider(1:300, show_value=true))"
# ╔═╡ 372e190b-0344-4001-9c72-171610b41e9c
img_noisy = abs.(simulate(img, p, r2));
# ╔═╡ b17cd903-bce6-4291-9d71-ea8a8422295a
gray_show(img_noisy, set_one=true)
# ╔═╡ fc91074a-99fc-4a13-802c-740872bf7502
md"
Mean is $(round(mean(img_noisy), sigdigits=3))
Maximum is $(round(maximum(img_noisy), sigdigits=3))
Minimum is $(round(minimum(img_noisy), sigdigits=3))
"
# ╔═╡ fac6c740-93c1-4e12-9610-f3c621f84043
md"## 5.1 Undersampling
Demonstration of aliasing.
Write a function `undersample` which only takes every `n-th` point of an array, keeping the first pixel of the array.
"
# ╔═╡ 15aa846c-7b81-4094-8b5b-dd6575e82afe
img3 = Float32.(testimage("resolution_test_512"));
# ╔═╡ be3b0811-c05a-4e78-9d34-8c52195dd7d0
"""
undersample(x, factor)
```julia-repl
julia> undersample([1,2,3,4,5,6,7,8,9,10], 2)
5-element Vector{Int64}:
1
3
5
7
9
julia> undersample([1,2,3,4,5,6,7,8,9,10], 3)
4-element Vector{Int64}:
1
4
7
10
julia> undersample([1,2,3,4,5,6,7,8,9,10], 4)
3-element Vector{Int64}:
1
5
9
```
"""
function undersample(x, factor)
return ones(size(x) .÷ factor)
end
# ╔═╡ e5188537-b955-4637-bab5-a76b6d74c189
PlutoTest.@test undersample([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 2) == [1, 3, 5, 7, 9]
# ╔═╡ c2b3cdc7-60b5-4287-8268-f95184d452b7
PlutoTest.@test undersample([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 3) == [1, 4, 7, 10]
# ╔═╡ 4e3a272e-0280-4e27-b014-3f9da5ac8039
PlutoTest.@test undersample([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 4) == [1, 5, 9]
# ╔═╡ a76e4081-7c52-4548-b513-ba03954378db
# ╔═╡ f5fbd0da-fe4c-4162-be5e-f9d9d5b27d87
md"## 5.1 Undersampling Visualization
In this plot you can see a nicely oversampled $\sin$ curve.
"
# ╔═╡ a07cb23e-6fca-4c42-a6b6-a525b2b5e232
xs2 = range(0, 10π, length=500)
# ╔═╡ c373046e-04b3-4159-80c0-8fbec141833f
f2(x) = sin(5 * x)
# ╔═╡ 9c1ed9aa-9d52-425d-892d-af3b8a981663
ys2 = f2.(xs2);
# ╔═╡ 75006409-0338-470d-a41d-fa1aba5c4fca
plot(xs2, ys2, xlabel="x pos", ylabel="amplitude")
# ╔═╡ f2b3d8a0-4bea-4443-9bc3-37258106c3d3
md"
Now we only take each `undersample_factor` data point of the series.
`undersample_factor` = $(@bind undersample_factor Slider(1:40, show_value=true))
Analyze what happens once you're below the Nyquist frequency.
See further what happens if you undersample then even more.
"
# ╔═╡ f6b06d6e-b268-45d6-a7e4-d8614cd0884d
y2_undersampled = undersample(ys2, undersample_factor);
# ╔═╡ eab1fe28-f336-4b33-89a3-4c65ed2a818f
plot(fftfreq(length(y2_undersampled))[1:length(y2_undersampled)÷2],abs.(rfft(y2_undersampled))[1:end-1], xlabel="Positive Frequencies", ylabel="Absolute Value of Frequency amplitude", mark="-*")
# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
FourierTools = "b18b359b-aebc-45ac-a139-9c0ccbb2871e"
ImageShow = "4e3cecfd-b093-5904-9786-8bbb286a6a31"
Noise = "81d43f40-5267-43b7-ae1c-8b967f377efa"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoTest = "cb4044da-4d16-4ffa-a6a3-8cad7f73ebdc"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
[compat]
Colors = "~0.12.8"
FFTW = "~1.5.0"
FourierTools = "~0.3.7"
ImageShow = "~0.3.6"
Noise = "~0.3.2"
Plots = "~1.36.2"
PlutoTest = "~0.2.2"
PlutoUI = "~0.7.48"
TestImages = "~1.7.1"
"""
# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised
julia_version = "1.8.2"
manifest_format = "2.0"
project_hash = "948be8f4bfb81fb5b88dd00248ba84f862ce6d3f"
[[deps.AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "69f7020bd72f069c219b5e8c236c1fa90d2cb409"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.2.1"
[[deps.AbstractNFFTs]]
deps = ["LinearAlgebra", "Printf"]
git-tree-sha1 = "292e21e99dedb8621c15f185b8fdb4260bb3c429"
uuid = "7f219486-4aa7-41d6-80a7-e08ef20ceed7"
version = "0.8.2"
[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"
[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "195c5505521008abea5aee4f96930717958eac6f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.4.0"
[[deps.ArgCheck]]
git-tree-sha1 = "a3a402a35a2f7e0b87828ccabbd5ebfbebe356b4"
uuid = "dce04be8-c92d-5529-be00-80e4d2c0e197"
version = "2.3.0"
[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
version = "1.1.1"
[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
[[deps.AxisArrays]]
deps = ["Dates", "IntervalSets", "IterTools", "RangeArrays"]
git-tree-sha1 = "1dd4d9f5beebac0c03446918741b1a03dc5e5788"
uuid = "39de3d68-74b9-583c-8d2d-e117c070f3a9"
version = "0.4.6"
[[deps.BangBang]]
deps = ["Compat", "ConstructionBase", "Future", "InitialValues", "LinearAlgebra", "Requires", "Setfield", "Tables", "ZygoteRules"]
git-tree-sha1 = "7fe6d92c4f281cf4ca6f2fba0ce7b299742da7ca"
uuid = "198e06fe-97b7-11e9-32a5-e1d131e6ad66"
version = "0.3.37"
[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
[[deps.Baselet]]
git-tree-sha1 = "aebf55e6d7795e02ca500a689d326ac979aaf89e"
uuid = "9718e550-a3fa-408a-8086-8db961cd8217"
version = "0.1.1"
[[deps.BasicInterpolators]]
deps = ["LinearAlgebra", "Memoize", "Random"]
git-tree-sha1 = "3f7be532673fc4a22825e7884e9e0e876236b12a"
uuid = "26cce99e-4866-4b6d-ab74-862489e035e0"
version = "0.7.1"
[[deps.BitFlags]]
git-tree-sha1 = "629c6e4a7be8f427d268cebef2a5e3de6c50d462"
uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35"
version = "0.1.6"
[[deps.Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"
[[deps.CEnum]]
git-tree-sha1 = "eb4cb44a499229b3b8426dcfb5dd85333951ff90"
uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
version = "0.4.2"
[[deps.Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.1+1"
[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "e7ff6cadf743c098e08fca25c91103ee4303c9bb"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.15.6"
[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "38f7a08f19d8810338d4f5085211c7dfa5d5bdd8"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.4"
[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"
[[deps.ColorSchemes]]
deps = ["ColorTypes", "ColorVectorSpace", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "1fd869cc3875b57347f7027521f561cf46d1fcd8"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.19.0"
[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "eb7f0f8307f71fac7c606984ea5fb2817275d6e4"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.4"
[[deps.ColorVectorSpace]]
deps = ["ColorTypes", "FixedPointNumbers", "LinearAlgebra", "SpecialFunctions", "Statistics", "TensorCore"]
git-tree-sha1 = "d08c20eef1f2cbc6e60fd3612ac4340b89fea322"
uuid = "c3611d14-8923-5661-9e6a-0046d554d3a4"
version = "0.9.9"
[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"
[[deps.Compat]]
deps = ["Dates", "LinearAlgebra", "UUIDs"]
git-tree-sha1 = "aaabba4ce1b7f8a9b34c015053d3b1edf60fa49c"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "4.4.0"
[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
version = "0.5.2+0"
[[deps.CompositionsBase]]
git-tree-sha1 = "455419f7e328a1a2493cabc6428d79e951349769"
uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b"
version = "0.1.1"
[[deps.ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "fb21ddd70a051d882a1686a5a550990bbe371a95"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.4.1"
[[deps.ContextVariablesX]]
deps = ["Compat", "Logging", "UUIDs"]
git-tree-sha1 = "25cc3803f1030ab855e383129dcd3dc294e322cc"
uuid = "6add18c4-b38d-439d-96f6-d6bc489c04c5"
version = "0.1.3"
[[deps.Contour]]
git-tree-sha1 = "d05d9e7b7aedff4e5b51a029dced05cfb6125781"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.6.2"
[[deps.DataAPI]]
git-tree-sha1 = "e08915633fcb3ea83bf9d6126292e5bc5c739922"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.13.0"
[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "d1fff3a548102f48987a52a2e0d114fa97d730f0"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.13"
[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"
[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
[[deps.DefineSingletons]]
git-tree-sha1 = "0fba8b706d0178b4dc7fd44a96a92382c9065c2c"
uuid = "244e2a9f-e319-4986-a169-4d1fe445cd52"
version = "0.1.2"
[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
[[deps.Distances]]
deps = ["LinearAlgebra", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "3258d0659f812acde79e8a74b11f17ac06d0ca04"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.7"
[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "c36550cb29cbe373e95b3f40486b9a4148f89ffd"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.9.2"
[[deps.Downloads]]
deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
version = "1.6.0"
[[deps.Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "bad72f730e9e91c08d9427d5e8db95478a3c323d"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.4.8+0"
[[deps.FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"
[[deps.FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Pkg", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "74faea50c1d007c85837327f6775bea60b5492dd"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.4.2+2"
[[deps.FFTW]]
deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
git-tree-sha1 = "90630efff0894f8142308e334473eba54c433549"
uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
version = "1.5.0"
[[deps.FFTW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c6033cc3892d0ef5bb9cd29b7f2f0331ea5184ea"
uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
version = "3.3.10+0"
[[deps.FLoops]]
deps = ["BangBang", "Compat", "FLoopsBase", "InitialValues", "JuliaVariables", "MLStyle", "Serialization", "Setfield", "Transducers"]
git-tree-sha1 = "ffb97765602e3cbe59a0589d237bf07f245a8576"
uuid = "cc61a311-1640-44b5-9fba-1b764f453329"
version = "0.2.1"
[[deps.FLoopsBase]]
deps = ["ContextVariablesX"]
git-tree-sha1 = "656f7a6859be8673bf1f35da5670246b923964f7"
uuid = "b9860ae5-e623-471e-878b-f6a53c775ea6"
version = "0.1.1"
[[deps.FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "7be5f99f7d15578798f338f5433b6c432ea8037b"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.16.0"
[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"
[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"
[[deps.Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.93+0"
[[deps.Formatting]]
deps = ["Printf"]