forked from fudanxu/SAR-GGCS
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathSea1.m
128 lines (123 loc) · 4.76 KB
/
Sea1.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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
clc
clear
close all
load('.\Data\sea1.mat');
hh = sea1;
c1 = -40; c2 = 10;
c11 = -5; c22 = 3;
N = 100;
% *************************
addpath(genpath(pwd));
% ********data processing
hh = double(hh);
sz = size(hh);
hh = hh(1:sz(1)-mod(sz(1),2),1:sz(2)-mod(sz(2),2));
[row,col] = size(hh);
zr = real(hh);
zi = imag(hh);
z = hh.*conj(hh);
za = sqrt(z);
%**************************parameter estimate*****************************************************
m = mean(mean(sqrt(hh.*conj(hh))));
sigma = 4*m^2/pi;
A = linspace(min(sqrt(hh(:).*conj(hh(:)))), max(sqrt(hh(:).*conj(hh(:)))), 100);
fprintf('The eatimated parameter is sigma = %.4f\n',sigma);
n = 1;
b = sqrt(sigma./2);
P_A = raylpdf(A,b);
counts = hist(za(:),A);
P_A_real = counts./numel(za)./(A(end)-A(end-1));
%***************estimate the correlation of real SAR image***************************
rho_sr = autocorr2d(zr);
rho_si = autocorr2d(zi);
rho_s = (rho_sr+rho_si)./2;
rho_sa = autocorr2d(za);
rho_sI = autocorr2d(za.^2);
% ********************compute the Convolution kernel of real/imag image********************************
rho_z_mean = rho_s;
rho_z = rho_z_mean;
Rho_z = ifftshift(fft2(fftshift(rho_z)));
Rho_z(imag(Rho_z)<1e-10) = real(Rho_z(imag(Rho_z)<1e-10));
Hr = sqrt(abs(Rho_z));
hr = ifftshift(ifft2(fftshift(Hr)));
hr = Select_center_win(hr,20);
% ********************compute the Convolution kernel of Amplitude image ********************************
rho_za = rho_sa;
Rho_A = ifftshift(fft2(fftshift(rho_za)));
Rho_A(imag(Rho_A)<1e-10) = real(Rho_A(imag(Rho_A)<1e-10));
Ha = sqrt(Rho_A);
ha = ifftshift(ifft2(fftshift(Ha)));
ha(imag(ha)<1e-10) = real(ha(imag(ha)<1e-10));
% ********************compute the Convolution kernel of Intensity image ********************************
rho_zI = rho_sI;
Rho_I = ifftshift(fft2(fftshift(rho_zI)));
Rho_I(imag(Rho_I)<1e-10) = real(Rho_I(imag(Rho_I)<1e-10));
HI = sqrt(Rho_I);
hI = ifftshift(ifft2(fftshift(HI)));
hI(imag(hI)<1e-10) = real(hI(imag(hI)<1e-10));
%**************single-point distribution*******************************
%ÈðÀû·Ö²¼pdf
b = sqrt(sigma./2);
P_A = raylpdf(A,b);
%***********Coherent scatterer model*********************************
s = sqrt(sigma./2);
%*******************************Generate scatterer number N******************************************
%****************************************Generate the truth Scatterer number***************************************
sigmai = s./sqrt(N);
% N_truth = za.^2./(2*sigmai^2);
% N_truth = ceil(mean(N_truth(:)));
%
% fprintf('The fitted parameters of N_truth is N_truth = %.4f\n',N_truth);
% fprintf('The estimated parameters of N is N = %.4f\n',N);
%nerate the gaussian scattering field
rng(1000000);
R_simu = normrnd(0,s./sqrt(N),row,col,N);
I_simu = normrnd(0,s./sqrt(N),row,col,N);
% ********************generate the correlated gaussian scattering field********************************
h0 = waitbar(0,'Please wait...');
for k = 1:max(max(N))
R_simu(:,:,k) = imfilter(R_simu(:,:,k),hr,'replicate');
[m,sigma] = normfit(reshape(R_simu(:,:,k),1,[]));
R_simu(:,:,k) = (R_simu(:,:,k)-m)./sigma;
R_simu(:,:,k) = R_simu(:,:,k).*s./sqrt(N);
I_simu(:,:,k) = imfilter(I_simu(:,:,k),hr,'replicate');
[m,sigma] = normfit(reshape(I_simu(:,:,k),1,[]));
I_simu(:,:,k) = (I_simu(:,:,k)-m)./sigma;
I_simu(:,:,k) = I_simu(:,:,k).*s./sqrt(N);
waitbar(k/max(max(N)),h0);
end
close(h0)
%***obtain the real image and imaginary image ************
R_simu_sum = zeros(row,col);
I_simu_sum = zeros(row,col);
% *********first N summation**********
for k1 = 1:row
for k2 = 1:col
R_simu_sum(k1,k2,:) = sum(R_simu(k1,k2,1:N));
I_simu_sum(k1,k2,:) = sum(I_simu(k1,k2,1:N));
end
end
%***obtain Amplitude image, Intensity image ************
A_simu = sqrt(R_simu_sum.^2+I_simu_sum.^2);
theta = atan2(I_simu_sum,R_simu_sum);
II_simu = A_simu.^2;
counts = hist(za(:),A);
P_A_real = counts./numel(za)./(A(end)-A(end-1));
[counts,centers] = hist(atan2(zi(:),zr(:)),100);
P_Theta_real = counts./numel(zi)./(centers(end)-centers(end-1));
counts = hist(A_simu(:),A);
P_A_simu = counts./numel(A_simu)./(A(end)-A(end-1));
[counts,centers] = hist(theta(:),100);
P_Theta_simu = counts./numel(theta)./(centers(end)-centers(end-1));
A_corr = autocorr2d(A_simu);
rho_sa = autocorr2d(za);
A_corr_w = Select_center_win(A_corr,10);
rho_sa_w = Select_center_win(rho_sa,10);
%% £¨a£©real data
figure; imagesc(20.*log10(za(1:row,1:col))); axis equal tight off; colorbar; colormap('gray'); caxis([-40,10]);
tightfig;
print('-dtiff','-r300',['Results\','Sea1_real'])
%% £¨b£©simulated data
figure; imagesc(20.*log10(A_simu)); axis equal tight off; colorbar; colormap('gray'); caxis([-40,10]);
tightfig;
print('-dtiff','-r300',['Results\','Sea1_simulated'])