囚徒_风云云检测算法改进

function mask = code(ref_b2,ref_b3,ref_b4,ref_b5,tmp_7,tmp_9,tmp_13,tmp_15,SC,height,mask_lan)
%算法实现
% 此处提供详细说明
sz=size(ref_b2);
temp=ref_b4*0;
temp(temp~=-999.0)=1;
raio=ref_b3 ./ref_b2;
%可信矩阵准备
mat_15=T_mat(tmp_15,224,228,"lt");
mat_9=T_mat(tmp_9,215,225,"lt");
mat_2=T_mat(ref_b2,0.38,0.34,"gt");
mat_3=T_mat(ref_b3,0.23,0.1,"gt");
mat_13=T_mat(tmp_13,280,290,"lt");
mat_4=zeros(sz(1),sz(2));
for i=1:sz(1)
for j=1:sz(2)
if (ref_b4(i,j) >0.06)&&(height(i,j)<3000)&&(temp(i,j)==1)
mat_4(i,j)=0;
elseif (ref_b4(i,j) <0.04)&&(temp(i,j)==1) || ((ref_b4(i,j) >0.06)&&(height(i,j)>3000))
mat_4(i,j)=1;
else
mat_4(i,j)=(0.06-ref_b4(i,j))/0.02;
end
end
end
mat_ra=zeros(sz(1),sz(2));
for i=1:sz(1)
for j=1:sz(2)
if abs(SC(i,j)) < 20
if (raio(i,j) >0.92)
mat_ra(i,j)=0;
elseif raio(i,j) <0.9
mat_ra(i,j)=1;
else
mat_ra(i,j)=(0.92-raio(i,j))/(0.02);
end
else
if raio(i,j) > 0.8
mat_ra(i,j)=0;
elseif raio(i,j) <0.78
mat_ra(i,j)=1;
else
mat_ra(i,j)=(0.8-raio(i,j))/(0.02);
end
end
end
end
tmp_7_13=zeros(sz(1),sz(2));
tmp_7_13=tmp_7-tmp_13;
mat_7_13=zeros(sz(1),sz(2));
for i=1:sz(1)
for j=1:sz(2)
if abs(SC(i,j)) < 20
if tmp_7_13(i,j) > 35.5
mat_7_13(i,j)=0;
elseif tmp_7_13(i,j) <33
mat_7_13(i,j)=1;
else
mat_7_13(i,j)=(35.5-tmp_7_13(i,j))/(2.5);
end
else
if tmp_7_13(i,j) > 19
mat_7_13(i,j)=0;
elseif tmp_7_13(i,j) <17
mat_7_13(i,j)=1;
else
mat_7_13(i,j)=(19-tmp_7_13(i,j))/(2.0);
end
end
end
end
va=zeros(sz(1),sz(2));
va=(ref_b2-ref_b5)./(ref_b2+ref_b5);
%沙漠处理
for i=1:sz(1)
for j=1:sz(2)
if (va(i,j) < -0.04)&&(mask_lan(i,j) == 1)
mat_7_13(i,j)=1;
end
end
end
%陆地
%高云
p_15_lan=zeros(sz(1),sz(2));
p_9_lan=zeros(sz(1),sz(2));
p_4_lan=zeros(sz(1),sz(2));
%高云可信度
F1_lan=zeros(sz(1),sz(2));
temp1_lan=zeros(1,3);
%中云
p_2_lan=zeros(sz(1),sz(2));
%可信度
F2_lan=zeros(sz(1),sz(2));
%低云
p_7_13_lan=zeros(sz(1),sz(2));
F3_lan=zeros(sz(1),sz(2));
%%海洋
%高云
p_15_sea=zeros(sz(1),sz(2));
p_9_sea=zeros(sz(1),sz(2));
p_4_sea=zeros(sz(1),sz(2));
p_13_sea=zeros(sz(1),sz(2));
%高云可信度
F1_sea=zeros(sz(1),sz(2));
temp1_sea=zeros(1,4);
%中云
p_3_sea=zeros(sz(1),sz(2));
p_ra_sea=zeros(sz(1),sz(2));
%可信度
F2_sea=zeros(sz(1),sz(2));
temp2_sea=zeros(1,2);
%低云
p_7_13_sea=zeros(sz(1),sz(2));
F3_sea=zeros(sz(1),sz(2));
%%
W=2;
for i=1+W:sz(1)-W
for j=1+W:sz(2)-W
if mask_lan (i,j)==1
%高云监测
p_15_lan(i,j)=mat_15(i,j);
if (tmp_15(i,j)>224)&&(tmp_15(i,j)<228)
value=inter(mat_15,i,j,mat_15(i,j),1);
p_15_lan(i,j)=value;
end
p_9_lan(i,j)=mat_9(i,j);
if (tmp_9(i,j)>215)&&(tmp_9(i,j)<225)
value=inter(mat_9,i,j,mat_9(i,j),1);
p_9_lan(i,j)=value;
end
p_4_lan(i,j)=mat_4(i,j);
if (ref_b4(i,j)>0.04)&&(ref_b4(i,j)<0.06)
value=inter(mat_4,i,j,mat_4(i,j),1);
p_4_lan(i,j)=value;
end
%高云的可信度
temp1_lan(3)=p_15_lan(i,j);
temp1_lan(1)=p_9_lan(i,j);
temp1_lan(2)=p_4_lan(i,j);
F1_lan(i,j)=min(temp1_lan);
%中云
p_2_lan(i,j)=mat_2(i,j);
if (ref_b2(i,j)>0.34)&&(ref_b2(i,j)<0.38)
value=inter(mat_2,i,j,mat_2(i,j),1);
p_2_lan(i,j)=value;
end
%可信度
F2_lan(i,j)=p_2_lan(i,j);
%低云
p_7_13_lan(i,j)=mat_7_13(i,j);
if (tmp_7_13(i,j)>17)&&(tmp_7_13(i,j)<19)
value=inter(mat_7_13,i,j,mat_7_13(i,j),1);
p_7_13_lan(i,j)=value;
end
%可信度
F3_lan(i,j)=p_7_13_lan(i,j);
else %海
%高云监测
p_15_sea(i,j)=mat_15(i,j);
if (tmp_15(i,j)>224)&&(tmp_15(i,j)<228)
value=inter(mat_15,i,j,mat_15(i,j),1);
p_15_sea(i,j)=value;
end
p_9_sea(i,j)=mat_9(i,j);
if (tmp_9(i,j)>215)&&(tmp_9(i,j)<225)
value=inter(mat_9,i,j,mat_9(i,j),1);
p_9_sea(i,j)=value;
end
p_4_sea(i,j)=mat_4(i,j);
if (ref_b4(i,j)>0.04)&&(ref_b4(i,j)<0.06)
value=inter(mat_4,i,j,mat_4(i,j),1);
p_4_sea(i,j)=value;
end
p_13_sea(i,j)=mat_13(i,j);
if (tmp_13(i,j)>280)&&(tmp_13(i,j)<290)
value=inter(mat_7_13,i,j,mat_7_13(i,j),1);
p_13_sea(i,j)=value;
end
%高云的可信度
temp1_sea(4)=p_15_sea(i,j);
temp1_sea(1)=p_9_sea(i,j);
temp1_sea(2)=p_4_sea(i,j);
temp1_sea(3)=p_13_sea(i,j);
F1_sea(i,j)=min(temp1_sea);
%中云
p_3_sea(i,j)=mat_3(i,j);
if (ref_b3(i,j)>0.1)&&(ref_b3(i,j)<0.23)
value=inter(mat_3,i,j,mat_3(i,j),1);
p_3_sea(i,j)=value;
end
p_ra_sea(i,j)=mat_ra(i,j);
if abs(SC(i,j))<20
if (raio(i,j)>0.9)&&(raio(i,j)<0.92)
value=inter(mat_ra,i,j,mat_ra(i,j),1);
p_ra_sea(i,j)=value;
end
else
if (raio(i,j)>0.78)&&(raio(i,j)<0.8)
value=inter(mat_ra,i,j,mat_ra(i,j),1);
p_ra_sea(i,j)=value;
end
end
%可信度
temp2_sea(2)=p_3_sea(i,j);
temp2_sea(1)=p_ra_sea(i,j);
F2_sea(i,j)=min(temp2_sea);
%低云
p_7_13_sea(i,j)=mat_7_13(i,j);
if abs(SC(i,j))<20
if (tmp_7_13(i,j)>33)&&(tmp_7_13(i,j)<35.5)
value=inter(mat_7_13,i,j,mat_7_13(i,j),1);
p_7_13_sea(i,j)=value;
end
else
if (raio(i,j)>17)&&(raio(i,j)<19)
value=inter(mat_7_13,i,j,mat_7_13(i,j),1);
p_7_13_sea(i,j)=value;
end
end
%可信度
F3_sea(i,j)=p_7_13_sea(i,j);
end
end
end
%可信度海洋F_sea
F_sea=zeros(sz(1),sz(2));
F_sea=nthroot((F1_sea.*F2_sea.*F3_sea),3);
%针对耀斑区的晴空优化
for i=1:sz(1)
for j=1:sz(2)
if abs(SC(i,j))<20
if F_sea(i,j)>0.85
F_sea(i,j)=1;
end
end
end
end
%可信度陆地F_lan
F_lan=zeros(sz(1),sz(2));
F_lan=nthroot((F1_lan.*F2_lan.*F3_lan),3);
%可信度F
F=zeros(sz(1),sz(2));
F=F_lan+F_sea;
%晴空修复
ndvi=(ref_b3-ref_b2)./(ref_b3+ref_b2);
for i=1:sz(1)
for j=1:sz(2)
if (ndvi(i,j) > 0.4)||(ndvi(i,j) <-0.18)
F(i,j)=1;
end
end
end
mask=mask_co(ref_b2,sz,F,mask_lan);
disp("code OKAY!")
end