-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmyIterativeDenoising.m
executable file
·113 lines (96 loc) · 3.83 KB
/
myIterativeDenoising.m
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
%% Making use of residuals in image denoising
% We aim to demonstrate the use of residuals in the evaluation and
% improvement of the various denoising algorithms we have studied in this
% course.
% The following are the meanings of the symbols used:
% 1. $\mathbf{x}$: Original image
% 2. $\mathbf{n}$: Additive white Gaussian noise
% 3. $\sigma^2$: Noise power
% 4. $\mathbf{y}$: Noisy image
% 5. $\mathbf{d}$: Denoised image
% 6. $\mathbf{r}$: Residual of the image $(\mathbf{d}-\mathbf{y})$
%%
% We prepare the field by loading the data and adding noise to it.
x_img = imread("../images/lena_256.tif");
x_img = double(x_img);
sigma2 = 400;
y_img = x_img + sqrt(sigma2) * randn(size(x_img));
figure(1);
imshow(mat2gray(x_img));
figure(2);
imshow(mat2gray(y_img));
%% Applying first iteration of denoising
% Here we choose the algorithm being used to denoise the image $y$.
d_img = myBilateralFilter(y_img, 6, 15, 20);
figure(3);
imshow(mat2gray(d_img));
%%
% Now we calculate the residual image.
r_img = y_img - d_img;
%%
% Next we find out the no-reference peak signal-to-noise ratio estimate
% using the estimate provided in the paper. The estimate is given by:
%
% $$\mathrm{PSNR}(\mathbf{x}, \mathbf{d}) = 10 \log_{10}\left(
% \dfrac{L^2}{\hat{\mathrm{MSE}}(\mathbf{x}, \mathbf{d})} \right)$$
%
% The mean-squared error estimate between the denoised image and the
% original image is given by:
%
% $$ \hat{\mathrm{MSE}}(\mathbf{x}, \mathbf{d})} = \overbar{\mathbf{r^2}} +
% \sigma^2 - 2\,\min\left( \overbar{\mathbf{r^2}}, s_{\mathbf{yr}},
% \sigma^2} \right) $$
%
% where $\overbar{\mathbf{r^2}}$ is the sample mean of $\mathbf{r^2}$ and
% $s_{\mathbf{yr}}$ is the sample covariance of $\mathbf{y}$ and
% $\mathbf{r}$.
r2_bar = mean(r_img.^2, 'all');
s_yr = mean((y_img - mean(y_img, 'all')).*(r_img - mean(r_img, 'all')), 'all');
mse = r2_bar + sigma2 - 2 * min(r2_bar, min(s_yr, sigma2));
psnr_val = 10*log10((255)^2/mse);
%% Iterative Denoising Algorithm
% We use an iterative denoising algorithm to denoise our image. The result
% is guaranteed to be atleast as good as the initially denoised image.
num_iters = 5;
y_img_array = zeros([size(x_img) num_iters]);
d_img_array = zeros([size(x_img) num_iters]);
r_img_array = d_img_array;
rd_img_array = d_img_array;
y_img_array(:,:,1) = y_img;
y_quality = zeros(num_iters, 1);
y_quality(1) = quality(y_img, d_img, sigma2);
d_quality = zeros(num_iters, 1);
for iter = 2:num_iters
d_img_array(:,:,iter) = myBilateralFilter(y_img_array(:,:,iter-1), 6, 15, 20);
d_quality(iter) = quality(y_img, d_img_array(:,:,iter), sigma2);
r_img_array(:,:,iter) = y_img - d_img_array(:,:,iter);
rd_img_array(:,:,iter) = myPCADenoising(r_img_array(:,:,iter));
y_img_array(:,:,iter) = d_img_array(:,:,iter) + rd_img_array(:,:,iter);
y_quality(iter) = quality(y_img, y_img_array(:,:,iter), sigma2);
figure(8*iter+1);
imshow(mat2gray(d_img_array(:,:,iter)));
[pearsonsmatrix] = pearsonscoeff(d_img_array(:,:,iter),y_img-d_img_array(:,:,iter));
figure(8*iter+2),
imshow(abs(pearsonsmatrix), [0,1.0]);
figure(8*iter+4),
imshow(y_img-d_img_array(:,:,iter), [-100,100]);
[kstestimgthresh] = kstestfun((y_img-d_img_array(:,:,iter))/20);
[autocorrimg] = autocorr(y_img-d_img_array(:,:,iter));
figure(8*iter+6),
imshow(kstestimgthresh, [0,1.0]);
figure(8*iter+7),
imshow(autocorrimg, [0,49.0]);
end
y_quality
d_quality
[y_qual_best, y_iter_max] = max(y_quality(2:end));
y_iter_max = y_iter_max + 1;
[d_qual_best, d_iter_max] = max(d_quality);
d_best = d_img_array(:,:,d_iter_max);
if (y_qual_best > d_qual_best)
d_best = y_img_array(:,:,y_iter_max);
d_qual_best = y_qual_best;
end
figure(4);
imshow(mat2gray(d_best));
fprintf('Final quality value = %.4f\n\n', d_qual_best);