%-----子函数adp_spa_denoise_guassian---------
%-------求出整个差值矩阵，再求这个矩阵得均值和方差相加后得到域值M，再来做噪声判断-------------
function [x,noise]=adp_spa_denoise_guassian(t,r)
[dep,wide]=size(t);
x=ones(size(t));
noise=zeros(size(t));
delta=ones(size(t));
for i=3:dep-2
for j=3:wide-2
b=sort([abs(t(i-2,j-2)-t(i,j)) abs(t(i-2,j-1)-t(i,j)) abs(t(i-2,j)-t(i,j)) abs(t(i-2,j+1)-t(i,j)) abs(t(i-2,j+2)-t(i,j)) abs(t(i-1,j-2)-t(i,j)) abs(t(i-1,j-1)-t(i,j)) abs(t(i-1,j)-t(i,j)) abs(t(i-1,j+1)-t(i,j)) abs(t(i-1,j+2)-t(i,j)) abs(t(i,j-2)-t(i,j)) abs(t(i,j-1)-t(i,j)) abs(t(i,j)-t(i,j)) abs(t(i,j+1)-t(i,j)) abs(t(i,j+2)-t(i,j)) abs(t(i+1,j-2)-t(i,j)) abs(t(i+1,j-1)-t(i,j)) abs(t(i+1,j)-t(i,j)) abs(t(i+1,j+1)-t(i,j)) abs(t(i+1,j+2)-t(i,j)) abs(t(i+2,j-2)-t(i,j)) abs(t(i+2,j-1)-t(i,j)) abs(t(i+2,j)-t(i,j)) abs(t(i+2,j+1)-t(i,j)) abs(t(i+2,j+2)-t(i,j))]);
delta(i,j)=(b(4)+b(5)+b(6)+b(7))/4;
end
end
for i=3:dep-2
delta(i,1)=delta(i,3);
delta(i,2)=delta(i,3);
delta(i,wide)=delta(i,wide-2);
delta(i,wide-1)=delta(i,wide-2);
end
delta(1,:)=delta(3,:);
delta(2,:)=delta(3,:);
delta(dep,:)=delta(dep-2,:);
delta(dep-1,:)=delta(dep-2,:);
delta_mean=0;
N=wide*dep;
for i=1:dep
for j=1:wide
delta_mean=delta_mean+delta(i,j);
end
end
delta_mean=delta_mean/N;
delta_var=0;
for i=1:dep
for j=1:wide
delta_var=delta_var+(delta(i,j)-delta_mean);
end
end
delta_var=delta_var/sqrt(N);
M=(delta_mean+delta_var)*r;
for i=3:dep-2
for j=3:wide-2
cnt=0;
for m=-2:2
for n=-2:2
if(delta(i,j)>M)
cnt=cnt+1;
end
end
end
if(cnt>=16)
x(i,j)=mean([t(i-2,j-2) t(i-2,j-1) t(i-2,j) t(i-2,j+1) t(i-2,j+2) t(i-1,j-2) t(i-1,j-1) t(i-1,j) t(i-1,j+1) t(i-1,j+2) t(i,j-2) t(i,j-1) t(i,j) t(i,j+1) t(i,j+2) t(i+1,j-2) t(i+1,j-1) t(i+1,j) t(i+1,j+1) t(i+1,j+2) t(i+2,j-2) t(i+2,j-1) t(i+2,j) t(i+2,j+1) t(i+2,j+2)]);
noise(i,j)=1;
else
x(i,j)=0.5*t(i,j)+0.5*(t(i-1,j-1)+t(i-1,j)+t(i-1,j+1)+t(i,j-1)+t(i,j+1)+t(i+1,j-1)+t(i+1,j)+t(i+1,j+1))/8.0;
noise(i,j)=0;
end
end
end
for i=3:dep-2
x(i,1)=x(i,3);
x(i,2)=x(i,3);
x(i,wide-1)=x(i,wide-2);
x(i,wide)=x(i,wide-2);
end
x(1,:)=x(3,:);
x(2,:)=x(3,:);
x(dep-1,:)=x(dep-2,:);
x(dep,:)=x(dep-2,:);
