function check_matrix= ICE_MASK(ValidationPath,lat_1km,lon_1km,time,lat_resized,lon_resized)
%
%  check_matrix= ICE_MASK(ValidationPath,lat_1km,lon_1km,time,lon_resized,lat_resized)
% 
% functions to be used to extract the ice mask of SAR image from IMS data.
% As slowly as searching at IMS pixel-to pixel, this function seperates SAR
% image to several subimages and extracts the ice-covered matirx of every
% subimage. If the extracted matrix is not the same values, then search the
% ice-covered falg in the subimage pixel-to-pixel
%
% input:
%     ValidationPath: the path of IMS files
%     lat_1km: the latitude matrix of IMS files
%     lon_1km: the longitude matrix of IMS files
%     time: the sensing time of SAR image
%     lat_resized: the latitude matrix of SAR image
%     lon_resized: the longitude matrix of SAR image
% output:
%     check_matrix: the ice and land mask matrix

%% read the IMS data
time_str=num2str(time);
file0_temp = dir([ValidationPath,'\*',time_str,'*.tif']);   
file_temp_name=file0_temp.name;
tif_temp = [ValidationPath,'\',file_temp_name];         % get the file path
IceExtentDaily=double(imread(tif_temp));                % read IMS image

%% extract the subimage of IMS according to the four corner points of SAR image
lon_resized(lon_resized(:)>180)=lon_resized(lon_resized(:)>180)-360;       % keep the longitude range grom -180 to 180
lon_check=[lon_resized(1,1),lon_resized(1,end),lon_resized(end,1),lon_resized(end,end)];
% convert the longitude range from 0 to 360 if the image crosses the 0 degree longitude
if sum(lon_check>0)~=0||sum(lon_check>0)~=4
    lon_1km(lon_1km<0)= lon_1km(lon_1km<0)+360;
    lon_resized(lon_resized(:)<0)=lon_resized(lon_resized(:)<0)+360;
end

% extract latitude and longitude of four corner points of SAR image
lon_left_up=lon_resized(1,1);lon_left_down=lon_resized(end,1);lon_right_up=lon_resized(1,end);lon_right_down=lon_resized(end,end);
lon_4=[lon_left_up,lon_left_down,lon_right_up,lon_right_down];
lat_left_up=lat_resized(1,1);lat_left_down=lat_resized(end,1);lat_right_up=lat_resized(1,end);lat_right_down=lat_resized(end,end);
lat_4=[lat_left_up,lat_left_down,lat_right_up,lat_right_down];

% find the position of four points above in IMS 
[ lon_1_x1,lon_1_y1 ] = FindCornerLocation_rever( lon_1km,lat_1km,lon_4(1),lat_4(1));
[ lon_2_x1,lon_2_y1 ] = FindCornerLocation_rever( lon_1km,lat_1km,lon_4(2),lat_4(2));
[ lat_3_x1,lat_3_y1 ] = FindCornerLocation_rever( lon_1km,lat_1km,lon_4(3),lat_4(3));
[ lat_4_x1,lat_4_y1 ] = FindCornerLocation_rever( lon_1km,lat_1km,lon_4(4),lat_4(4));

x_matrix=[lon_1_x1,lon_2_x1,lat_3_x1,lat_4_x1];
y_matrix=[lon_1_y1,lon_2_y1,lat_3_y1,lat_4_y1];
x_min=min(x_matrix);x_max=max(x_matrix);
y_min=min(y_matrix);y_max=max(y_matrix);

% extract the image of IMS according to the locations of four corner points
IceExtentDaily_cut=IceExtentDaily(x_min-10:x_max+10,y_min-10:y_max+10);
imslon_1km_cut=lon_1km(x_min-10:x_max+10,y_min-10:y_max+10);
imslat_1km_cut=lat_1km(x_min-10:x_max+10,y_min-10:y_max+10);

%% match IMS image and SAR image pixel-to-pixel

[a,b]=size(lon_resized);
subsize=40;                     % the size of subimage 
check_matrix=zeros(a,b);
M=floor(a/subsize);             % the number of subimage in range direction
N=floor(b/subsize);             % the number of subimage in azimuth direction

for i=1:M:a
    for j=1:N:b
        % at the border of SAR image, the size of subimage will not be the
        % subsize set above, so the rows and cols need to be adjusted
        if i+M>a   mm=a-i;
        else mm=M;
        end
        if  j+N>b  nn=b-j;
        else nn=N;
        end  
        
        % extract the sub-subimage of IMS according to each SAR subimage
        lon_cut_s1=lon_resized(i:i+mm,j:j+nn);
        lat_cut_s1=lat_resized(i:i+mm,j:j+nn);
        lat_down=min(lat_cut_s1(:));lat_up=max(lat_cut_s1(:));
        lon_down=min(lon_cut_s1(:));lon_up=max(lon_cut_s1(:));
        
        [lon_(1),lat_(1)]=FindCornerLocation_rever( imslon_1km_cut,imslat_1km_cut,lon_down,lat_down);
        [lon_(2),lat_(2)]=FindCornerLocation_rever( imslon_1km_cut,imslat_1km_cut,lon_down,lat_up);
        [lon_(3),lat_(3)]=FindCornerLocation_rever( imslon_1km_cut,imslat_1km_cut,lon_up,lat_down);
        [lon_(4),lat_(4)]=FindCornerLocation_rever( imslon_1km_cut,imslat_1km_cut,lon_up,lat_up);
        lon_num=[min(lon_):max(lon_)];
        lat_num=[min(lat_):max(lat_)];
        
        IceExtentDaily_cut_cut=IceExtentDaily_cut(lon_num,lat_num);
        max_ICE=max(IceExtentDaily_cut_cut(:));
        min_ICE=min(IceExtentDaily_cut_cut(:));
        
        % judge whether the values of IMS sub-subimage are all the same
        if max_ICE==min_ICE
            switch(IceExtentDaily_cut_cut(1,1))
                case(1)             % all pixels of this subimage is open water
                    check_matrix(i:i+mm,j:j+nn)=1;
                case{2,4}           % all pixels of this subimage is is land
                    check_matrix(i:i+mm,j:j+nn)=-999;
                otherwise           % all pixels of this subimage is is sea ice
                    check_matrix(i:i+mm,j:j+nn)=-1;
            end
        else
            % the pixels of subimage are not the same, then match the subimage
            % with IceExtentDaily_cut_cut pixel to pixel
            lon_cut=imslon_1km_cut(lon_num,lat_num);
            lat_cut=imslat_1km_cut(lon_num,lat_num);
            lon_temp=lon_resized(i:i+mm,j:j+nn);
            lat_temp=lat_resized(i:i+mm,j:j+nn);
            check_matrix_temp=zeros(mm,nn);
            for m=1:mm
                for n=1:nn   
                    % find the position in IMS sub-subimage of every SAR pixels
                    [ x,y ] = FindCornerLocation_rever(lon_cut,lat_cut,lon_temp(m,n),lat_temp(m,n));
                    switch(IceExtentDaily_cut_cut(x,y))
                        case 1        % this pixel is open water
                            check_matrix_temp(m,n)=1;
                        case {2,4}    % this pixel is land
                            check_matrix_temp(m,n)=-999;
                        otherwise     % this pixel is sea ice
                            check_matrix_temp(m,n)=-1;  
                    end
                end
            end
            check_matrix(i:i+mm-1,j:j+nn-1)=check_matrix_temp;
        end
    end
end
end


