clc
close all
clear all;
%% 原始图像
f = imread('cameraman.tif');
f = double(f)/255;
figure,imshow(f);title('原始图像')
%% 降质观测图像1
l = f(1:30,1:30);
M = 256;
nhood = [5 5];
h1 = ones(nhood);
h1 = h1/25;
h1(M,M) = 0;
H1 = fft2(h1);
F = fft2(f);
G1 = F.*H1;
g1 = ifft2(G1);
sigma1 = (std2(l)^2/10^3)^(1/2);
g1 = imnoise(g1,'gaussian',0, sigma1^2);
figure,imshow(g1);title('低分辨图像')
%% 降质观测图像2
nhood = [7 7];
h2 = ones(nhood);
h2 = h2/49;
h2(M,M) = 0;
H2 = fft2(h2);
G2 = F.*H2;
g2 = ifft2(G2);
sigma2 = (std2(l)^2/10^4)^(1/2);
g2 = imnoise(g2,'gaussian',0, sigma2^2);
figure,imshow(g2);title('低分辨图像')
%% 降质观测图像3
nhood = [9 9];
h3 = ones(nhood);
h3 = h3/81;
h3(M,M) = 0;
H3 = fft2(h3);
G3 = F.*H3;
g3 = ifft2(G3);
sigma3 = (std2(l)^2/10^5)^(1/2);
g3 = imnoise(g3,'gaussian',0, sigma3^2);
figure,imshow(g3);title('低分辨图像')
%% 多通道LMMSE重构
g1 = real(g1);
g2 = real(g2);
g3 = real(g3);
G1 = G1.';
G2 = G2.';
G3 = G3.';
G1 = im2col(G1,[M M],'distinct');
G2 = im2col(G2,[M M],'distinct');
G3 = im2col(G3,[M M],'distinct');
G = [G1;G2;G3];
wlength = 3;
g1(M+wlength,M+wlength) = 0;
g2(M+wlength,M+wlength) = 0;
g3(M+wlength,M+wlength) = 0;
x = 0+eps*i
Diagonals = repmat(x,M^2,18);
for counter = 1:9
% Begin of formation of correlation matrices.
switch counter
case 1
for j = 0:M-1
for k = 0:M-1
sum = 0;
for n = 1:wlength
for m = 1:wlength
sum = sum+g1(n,m)*g1(n+j,m+k);
end;
end;
R (j+1,k+1) = 1/((wlength)^2)*sum;
end;
end;
case 2
for j = 0:M-1
for k = 0:M-1
sum = 0;
for n = 1:wlength
for m = 1:wlength
sum = sum+g2(n,m)*g1(n+j,m+k);
end;
end;
R (j+1,k+1) = 1/((wlength)^2)*sum;
end;
end;
case 3
for j = 0:M-1
for k = 0:M-1
sum = 0;
for n = 1:wlength
for m = 1:wlength
sum = sum+g3(n,m)*g1(n+j,m+k);
end;
end;
R (j+1,k+1) = 1/((wlength)^2)*sum;
end;
end;
case 4
for j = 0:M-1
for k = 0:M-1
sum = 0;
for n = 1:wlength
for m = 1:wlength
sum = sum+g1(n,m)*g2(n+j,m+k);
end;
end;
R (j+1,k+1) = 1/((wlength)^2)*sum;
end;
end;
case 5
for j = 0:M-1
for k = 0:M-1
sum = 0;
for n = 1:wlength
for m = 1:wlength
sum = sum+g2(n,m)*g2(n+j,m+k);
end;
end;
R (j+1,k+1) = 1/((wlength)^2)*sum;
end;
end;
case 6
for j = 0:M-1
for k = 0:M-1
sum = 0;
for n = 1:wlength
for m = 1:wlength
sum = sum+g3(n,m)*g2(n+j,m+k);
end;
end;
R (j+1,k+1) = 1/((wlength)^2)*sum;
end;
end;
case 7
for j = 0:M-1
for k = 0:M-1
sum = 0;
for n = 1:wlength
for m = 1:wlength
sum = sum+g1(n,m)*g3(n+j,m+k);
end;
end;
R (j+1,k+1) = 1/((wlength)^2)*sum;
end;
end;
case 8
for j = 0:M-1
for k = 0:M-1
sum = 0;
for n = 1:wlength
for m = 1:wlength
sum = sum+g2(n,m)*g3(n+j,m+k);
end;
end;
R (j+1,k+1) = 1/((wlength)^2)*sum;
end;
end;
case 9
for j = 0:M-1
for k = 0:M-1
sum = 0;
for n = 1:wlength
for m = 1:wlength
sum = sum+g3(n,m)*g3(n+j,m+k);
end;
end;
R (j+1,k+1) = 1/((wlength)^2)*sum;
end;
end;
end;% End of switch.
Pf = fft2(R);
Pf = Pf.';
column = im2col(Pf,[M M],'distinct');
Diagonals(:,counter) = column;
end;
for counter = 1:3
switch counter
case 1
H = H1;
case 2
H = H2;
case 3
H = H3;
end;
H = H.';
column = im2col(H,[M M],'distinct');
Diagonals(:,9+counter) = column;
Diagonals(:,12+counter) = conj(column);
end;
Diagonals(:,16) = sigma1^2;
Diagonals(:,17) = sigma2^2;
Diagonals(:,18) = sigma3^2;
I = linspace(1, M^2,M^2);
J = I;
A1 = sparse(I,J,Diagonals(:,1));
A2 = sparse(I,J+M^2,Diagonals(:,2));
A3 = sparse(I,J+2*M^2,Diagonals(:,3));
A4 = sparse(I+M^2,J,Diagonals(:,4));
A5 = sparse(I+M^2,J+M^2,Diagonals(:,5));
A6 = sparse(I+M^2,J+2*M^2,Diagonals(:,6));
A7 = sparse(I+2*M^2,J,Diagonals(:,7));
A8 = sparse(I+2*M^2,J+M^2,Diagonals(:,8));
A9 = sparse(I+2*M^2,J+2*M^2,Diagonals(:,9));
A1(3*M^2,3*M^2) = 0;
A2(3*M^2,3*M^2) = 0;
A3(3*M^2,3*M^2) = 0;
A4(3*M^2,3*M^2) = 0;
A5(3*M^2,3*M^2) = 0;
A6(3*M^2,3*M^2) = 0;
A7(3*M^2,3*M^2) = 0;
A8(3*M^2,3*M^2) = 0;
A = A1+A2+A3+A4+A5+A6+A7+A8+A9;
A1 = sparse(I,J,Diagonals(:,13));
A2 = sparse(I+M^2,J+M^2,Diagonals(:,14));
A3 = sparse(I+2*M^2,J+2*M^2,Diagonals(:,15));
A1(3*M^2,3*M^2) = 0;
A2(3*M^2,3*M^2) = 0;
B = A1+A2+A3;
column = Diagonals(:,10).* Diagonals(:,1);
column = column.* Diagonals(:,13);
column = column+ Diagonals(:,16);
Diagonals(:,1) = column;
column = Diagonals(:,10).* Diagonals(:,2);
column = column.* Diagonals(:,14);
Diagonals(:,2) = column;
column = Diagonals(:,10).* Diagonals(:,3);
column = column.* Diagonals(:,15);
Diagonals(:,3) = column;
column = Diagonals(:,11).* Diagonals(:,4);
column = column.* Diagonals(:,13);
Diagonals(:,4) = column;
column = Diagonals(:,11).* Diagonals(:,5);
column = column.* Diagonals(:,14);
column = column+ Diagonals(:,17);
Diagonals(:,5) = column;
column = Diagonals(:,11).* Diagonals(:,6);
column = column.* Diagonals(:,15);
Diagonals(:,6) = column;
column = Diagonals(:,12).* Diagonals(:,7);
column = column.* Diagonals(:,13);
Diagonals(:,7) = column;
column = Diagonals(:,12).* Diagonals(:,8);
column = column.* Diagonals(:,14);
Diagonals(:,8) = column;
column = Diagonals(:,12).* Diagonals(:,9);
column = column.* Diagonals(:,15);
column = column+ Diagonals(:,18);
Diagonals(:,9) = column;
A1 = sparse(I,J,Diagonals(:,1));
A2 = sparse(I,J+M^2,Diagonals(:,2));
A3 = sparse(I,J+2*M^2,Diagonals(:,3));
A4 = sparse(I+M^2,J,Diagonals(:,4));
A5 = sparse(I+M^2,J+M^2,Diagonals(:,5));
A6 = sparse(I+M^2,J+2*M^2,Diagonals(:,6));
A7 = sparse(I+2*M^2,J,Diagonals(:,7));
A8 = sparse(I+2*M^2,J+M^2,Diagonals(:,8)) ;
A9 = sparse(I+2*M^2,J+2*M^2,Diagonals(:,9));
A1(3*M^2,3*M^2) = 0;
A2(3*M^2,3*M^2) = 0;
A3(3*M^2,3*M^2) = 0;
A4(3*M^2,3*M^2) = 0;
A5(3*M^2,3*M^2) = 0;
A6(3*M^2,3*M^2) = 0;
A7(3*M^2,3*M^2) = 0;
A8(3*M^2,3*M^2) = 0;
C = A1+A2+A3+A4+A5+A6+A7+A8+A9;
C = inv(C);
A = A*B;
C = C*G;
F = A*C;
F1 = F(1:M^2,1);
F2 = F(M^2+1:2*M^2,1);
F3 = F(2*M^2+1:3*M^2,1);
F1 = col2im(F1,[M M],[M M],'distinct');
F1 = F1.';
F2 = col2im(F2,[M M],[M M],'distinct');
F2 = F2.';
F3 = col2im(F3,[M M],[M M],'distinct');
F3 = F3.';
ff1 = ifft2(F1);
ff2 = ifft2(F2);
ff3 = ifft2(F3);
ff1 = real(ff1);
ff2 = real(ff2);
ff3 = real(ff3);
%% 显示结果
figure,imshow(ff1)
figure;imshow(ff2)
figure;imshow(ff3)
