% This program is used to retrieve ocean significant wave height from the
% Sentinel-1 data in HH polarization using the developed BP neural network

% The program was created by: Ms. Ke Wu
% Affiliation: Aerospace Information Research Institution, Chinese Academy of Sciences
% email: wuke@radi.ac.cn
% Date of creation: 29nd, Oct., 2020

% clear
clc
%% User input
%input Sentinel-1 images
[Filename,File_path]=uigetfile('*.tiff','please choose HH-polarization SAR files');
vfilen=[File_path,Filename];
ValidationPath='K:\IMS\201911';                  % input the path of IMS files
% load('lon_1km.mat');							 % load the data of IMS longitude and latitude matrix
% load('lat_1km.mat');

if ~exist(vfilen)
    error('the file does not exist!');
end

%% read the Sentinel-1 images and its auxiliary files
Image=double(imread(vfilen));												% read raw Sentinel-1 images

idst=strfind(vfilen,'\');
S_filepath=vfilen(1:(max(idst)-12));
S_filename=vfilen((max(idst)+1):end-4);

xmln_HH=[S_filepath,'annotation\',S_filename,'xml'];                         % read Sentinel-1 annotation  files
[meta,~]= xml_read(xmln_HH);
xml_c_HH=[S_filepath,'annotation\calibration\calibration-',S_filename,'xml'];% read Sentinel-1 calibration files
[calibration_tree,~]=xml_read(xml_c_HH);
xml_n_HH=[S_filepath,'annotation\calibration\noise-',S_filename,'xml'];      % read Sentinel-1 noise files
[noise_tree,~]=xml_read(xml_n_HH);
clear Filename  File_path xmln_HH xml_c_HH xml_n_HH

%% radiation correction and geometric correction of Sentinel-1 image
time=str2num(S_filename(15:22));
[ sigma,lat2,lon2,inc2,noise_nesz ] = radiation_correction_and_removal_side( Image,calibration_tree,noise_tree,meta,time);
clear calibration_tree noise_tree Image

sigma_db=10*log10(sigma);                                                  % invert SAR image from linera unit to dB unit
noise_nesz(noise_nesz(:)<0)=0;
noise_db=10*log10(noise_nesz);                                             % calculate the NESZ
mul_db=sigma_db-noise_db;
mul_db(isinf(mul_db))=0;
clear noise_db noise_nesz sigma_db

% we consider SAR images at the north of 83 degree_north all covered
% by ice, so this program do not process these images
if min(lat2(:))>=83
    disp('lat>83');
    return
end

%% ice and land mask

Model_1024= ICE_MASK(ValidationPath,lat_1km,lon_1km,time,lat2,lon2);		% extract the ice mask matrix
mask_=Model_1024;
mask_(mask_(:)~=1)=0;
sigma_masked=sigma.*mask_;													% in this step, the pixels represented sea ice and land are all set to 0
clear mask_

% this program does not process the image with over 85% sea ice covered
[a,b]=size(Model_1024);
iceratio=sum(sum(Model_1024(:)>0))/(a*b);                                   % compute the ratio of ice covered in the whole image
if iceratio<0.15
   disp('ICE >0.85')
   return
end
clear a b iceratio

%%  extract parameters of swh retrieval
[rows,cols]=size(sigma_masked);
size_subscene=256;                                                          % size of subimage
dkx=2.0*pi./(size_subscene*16.0);
dky=2.0*pi./(size_subscene*20.0);

fft_sub=zeros(floor(rows/size_subscene*4),floor(cols/size_subscene*4),size_subscene/2,size_subscene/2);
sigma_sub1=zeros(floor(rows/size_subscene*4),floor(cols/size_subscene*4),size_subscene/2,size_subscene/2);
flag_sub=zeros(floor(rows/size_subscene*4),floor(cols/size_subscene*4));

% compute homogeneity value of each sub-subimage which size is subscene/2 and use fft to compute the image spectrum
m=0;
for i=size_subscene/4:size_subscene/4:rows-size_subscene/4
    m=m+1;n=0;
    for j=size_subscene/4:size_subscene/4:cols-size_subscene/4
        n=n+1;
        % extract SAR subimage
        sigma_subscene=sigma_masked(i-size_subscene/4+1:i+size_subscene/4,j-size_subscene/4+1:j+size_subscene/4);
        
        % remove the subimage of ice or land covered
        if sum(sum(sigma_subscene(:)==0|isnan(sigma_subscene(:))))==size_subscene.^2/4   
            continue
        end
        % calculate the homogeneity of subimage
        homo_sub=homogeneity(sigma_subscene);
        % remove the subimage which homogeneity value is over 1.05
        if homo_sub>1.05
            continue
        end
        clear homo_sub
        
        flag_sub(m,n)=1;                                                   % image_flag, flag=1 means the the subimage is homogeneous
        mean_subscene=nanmean(sigma_subscene(:));                          % compute the mean value of subimage
        intens=(sigma_subscene-mean_subscene)/mean_subscene;               % normalization of subimages
        intens(isnan(sigma_subscene(:)))=0;
        fft_subscene=ifftshift(abs(ifft2(intens)).^2)/dkx/dky;             % FFT
        fft_sub(m,n,:,:)=fft_subscene;
        sigma_sub1(m,n,:,:)=sigma_subscene;
        clear intens fft_subscene homo_sub sigma_subscene mean_subscene
    end
end
roww=m-2;  coll=n-2;
clear n m

% initialize the parameter matrix
sizee=size_subscene/4;                                                     % resampling spacing
lat_=zeros(roww,coll);  lon_=zeros(roww,coll);                             % the latitude and longitude matrix
sigma_submasked=zeros(roww,coll);                                          % the masked sigma matrix
sigma_unmasked=zeros(roww,coll);                                           % the unmasked sigma matrix                                                  %
cvar=zeros(roww,coll);                                                     % the cvar matrix
inc_256=zeros(roww,coll);                                                  % the incidence angle matrix
Mask=zeros(roww,coll);                                                     % the ice and land mask matrix
para_o=zeros(roww,coll,20);                                                % the spectrum parameters matrix (row*col*20)
flag=zeros(roww,coll);                                                     % the homogeneity flag matrix
flag_noise=zeros(roww,coll);                                               % the flag matrix of thermal noise
flag_unprocessed=zeros(roww,coll);                                         % the flag matrix of unprocessed pixel      

%% compute the parameters of the whole SAR images with the step of 2.56 km (64 pixels)
[h_o,index_1]=orth_func(256,256);                                          % compute 20 orthonormal functions
m=0;
for i=size_subscene/2:sizee:rows-size_subscene/2
    m=m+1;n=0;
    for j=size_subscene/2:sizee:cols-size_subscene/2
        n=n+1;
        
        lat_(m,n)=lat2(i,j);                                               % extract the center latitude of subimage
        lon_(m,n)=lon2(i,j);                                               % extract the center longitude of subimage
        Mask(m,n)=Model_1024(i,j);                                         % extract the mask flag of subimage
        sigma_unmasked(m,n)=nanmean(nanmean(sigma(i-sizee/2+1:i+sizee*1/2,...
            j-sizee/2+1:j+sizee*1/2)));                                    % extract the sigma of subimage
        inc_sub=inc2(i-size_subscene/2+1:i+size_subscene*1/2,...
            j-size_subscene/2+1:j+size_subscene*1/2);                      
        inc_256(m,n)=nanmean(inc_sub(:));                                  % extract the center incidence angle of subimage
        clear inc_sub
                
        % if the thermal noise is 3dB lower than mean sigma(dB) of subimage, the
        % flag value is set to 1
        flag_nesz=mul_db(i-sizee+1:i+sizee,j-sizee+1:j+sizee);
        if nanmean(flag_nesz(:))>3
            flag_noise(m,n)=1;
        end
        clear flag_nesz
        
        % if the subimage is all covered by sea ice or land, do not process it
        sigma_sub=sigma_masked(i-size_subscene/2+1:i+size_subscene/2,j-size_subscene/2+1:j+size_subscene/2);
        if sum(sigma_sub(:))==0   
            continue;
        end
        % compute the homgeneity value
        flag(m,n)=homogeneity(sigma_sub);
        clear sigma_sub
                
        % separate the subimage into four sub-subimages
        count=[m n;m+2 n;m n+2;m+2 n+2];
        count_sub=(count(:,2)-1)*floor(rows/size_subscene*4)+count(:,1);
        a=flag_sub(count_sub);
        count(a==0,:)=[];
        clear a count_sub
        % if more than 2 sub-subimages are imhomogeneous, the subimage will not be processed
        if size(count,1)<=2  
            flag_unprocessed(m,n)=1;
            continue
        end
        
        imfft=0;
        imsigma=[];
        % compute the average of sub-subimage and normalize it
        for kk=1:size(count,1)
            imsigma=[imsigma;squeeze(sigma_sub1(count(kk,1),count(kk,2),:,:))];
            imfft=imfft+squeeze(fft_sub(count(kk,1),count(kk,2),:,:));
        end
        fintens1=imfft/size(count,1);
        fintens1=imresize(fintens1,2);
        fintens1=fintens1/(sum(index_1(:))*dkx*dky);
        
        % compute 20 spectrum parameters from image spectrum
        for qq=1:20
            temp=fintens1.*h_o(:,:,qq);
            para_o(m,n,qq)=dkx*dky*sum(temp(isnan(temp(:))==0));
        end
        clear temp fintens1 imfft count
                
        % compute cvar
        imsigma(isnan(imsigma(:)))=0;
        mean_256=mean(imsigma(imsigma(:)~=0));
        mean_st=(imsigma-mean_256)/mean_256;                               % normalize the subimage
        stdm=std(mean_st(imsigma(:)>0));                                   % compute standard deviation 
        cvar(m,n)=stdm^2;
        sigma_submasked(m,n)=mean_256;
        clear imsigma stdm mean_st mean_256
    end
end
flag(flag(:)<=0)=0;
flag(flag_unprocessed==1)=0;

%% retrive swh
load('SWH_BP.mat');                                                  % load the trained BPNN model
m=roww;n=coll; clear roww coll

% convert the parameter matrixes to vectors
para_nrcs=reshape(sigma_submasked,m*n,1);
para_inc=cos(reshape(inc_256,m*n,1)*pi/180.0);
para_fft=reshape(para_o,m*n,20);
para_cvar=reshape(cvar,m*n,1);
input_swh=[para_nrcs para_inc para_cvar para_fft];

% compute swh using BPNN model and convert it to matrix
input_swhtest = mapminmax('apply',input_swh',inputps);
input_swhtest(isnan(input_swhtest(:)))=0;
outputswh=sim(net,input_swhtest);
swh_output = mapminmax('reverse',outputswh',outputps);
swh=reshape(swh_output,m,n);


%% draw swh map
ice=Mask;
ice(ice<0)=0;
SE=strel('square',3);
IMS_fcl1=imerode(ice,SE);
swh=swh.*IMS_fcl1;

% definite a new colormap
options.colormaps{1}.offset     = 0;
options.colormaps{1}.num        = 256;                                     % Number of colors
options.colormaps{1}.cmin       = -1    + options.colormaps{1}.offset;
options.colormaps{1}.cmax       = 255  + options.colormaps{1}.offset;
options.colormaps{1}.map        = [0.8 0.8 0.8;jet(256)];
cmap=options.colormaps{1}.map;

% compute the maximum and minimum of latitude and longitude
milo=min(lon_(:));malo=max(lon_(:));
mila=min(lat_(:));mala=max(lat_(:));
[a,b]=size(swh);

swh(isnan(swh(:)))=0;
h=figure('color','w');
set(h,'position',[520 110 750 600]);                           

% adjust the image position according to the length-width ratio of image
if b/a<=1.21
    set(gca,'units','normal','pos',[0.03 0.03 0.91 0.91],'PlotBoxAspectRatioMode','auto','DataAspectRatioMode','auto');
else
    set(gca,'units','normal','pos',[0.05 0.05 0.87 0.87],'PlotBoxAspectRatioMode','auto','DataAspectRatioMode','auto');
end

m_proj('Lambert Conformal Conic','long',[milo-0.08,malo+0.08],'lat',...
    [mila-0.08,mala+0.08]);                                                % map projection
m_grid('box','fancy','tickdir','in','Xaxislocation','top','fontsize',14);  % map grid
hold on 
pos_gca=get(gca,'position');
set(gca, 'LooseInset', [0,0,0,0.05]);                                      % make the image fill the figure
swh(swh(:)>=10|swh(:)<0|flag(:)>1.05|flag(:)==0|flag_noise(:)==0|ice(:)>=2)=nan;
swh(ice(:)<=0)=-1;                                                          % ice mask
h=m_pcolor(lon_,lat_,swh);                                                 % draw swh map
set(h,'edgecolor','none')
colormap(cmap);                                                            % set colormap
caxis([0,6]);                                                              % set the display range of swh
pos=[0.92,0.19,0.0213,pos_gca(4)*0.6667];                                  % set the position of colorbar
c=colorbar('FontSize',14,'position',pos);                                  % draw the colorbar
ylabel(c,'Significant Wave Height (m)','FontSize',16,'FontWeight',...
    'bold','Fontname','Times New Roman');    

% draw the data label
m_gshhs_h('patch',[0.86 0.66 0.41],'edgecolor','none');hold on             % draw the borderline
lon1=[lon_(1,1),lon_(1,end),lon_(end,end),lon_(end,1),lon_(1,1)];
lat1=[lat_(1,1),lat_(1,end),lat_(end,end),lat_(end,1),lat_(1,1)];
m_plot(lon1,lat1,'k-');hold on                                             % draw the four sides of image
m_grid('box','fancy','tickdir','in','Xaxislocation','top','fontsize',14);hold on 

% draw arrow of colorbar
x = pos(1,1); y = pos(1,2); w = pos(1,3); h = pos(1,4); 
ax3 = axes('Position',[x,y + h,w,0.05],'YAxisLocation','origin');
set(ax3,'ycolor','none','xcolor','none');
set(ax3,'XTick',[]);
fc = [1 2 3];
vt = [0.5 1;0 0;1 0;]; 
patch('Faces',fc,'Vertices',vt,'FaceColor',[0.5 0 0],'EdgeColor',[0.5 0 0]); hold on
ax4=axes('Position',[x-0.001,y-0.01,w+0.002,0.01+h/257*3],'YAxisLocation','origin');
set(ax4,'ycolor','none','xcolor','none');
set(ax4,'XTick',[]);hold on

%% print the image
file_name=S_filepath(end-72:end-1);
savename=[file_name(1:11),'_',file_name(18:32),'_',file_name(end-8:end-5)];

img_path=['SWH_',savename,'_V1.0.png'];
print(gcf, '-dpng', '-r600',img_path);