function [sigma0L_hh,lat2,lon2,inc2,noise_nesz]=radiation_correction_and_removal_side(Image,calibration_tree,noise_tree,meta,time)
%
%   [sigma0L_hh,lat2,lon2,inc2,noise_nesz]=radiation_correction_and_removal_side(Image,calibration_tree,noise_tree,meta,time)
%
% function of radiation correction and thermal noise removal of SAR image
% input:
%       image: the DN of SAR image
%       calibration_tree: the struct of calibration xml files
%       noise_tree: the struct of noise xml files
%       meta: the struct of annotation xml files
%       time: the sensing time of SAR image
% output: 
%       sigma0L_hh: the calibrated sigma of SAR
%       lat2: the latitude matrix of SAR
%       lon2: the longitude matrix of SAR
%       inc2: the incidence angle matrix of SAR
%       noise_nesc: equivalent noise of SAR image
%
[xf1,yf1]=size(Image);  
%% read and interpolate calibration vector
b2=1;
ll=length(calibration_tree.calibrationVectorList.calibrationVector);
co=calibration_tree.calibrationVectorList.calibrationVector(1).sigmaNought.ATTRIBUTE.count;
cali=zeros(ll,co);
for ii=1:ll    
    if ii==1
        cali(b2,1:co)=calibration_tree.calibrationVectorList.calibrationVector(ii).sigmaNought.CONTENT;        
        hh1=calibration_tree.calibrationVectorList.calibrationVector(ii).line;      
    else
        hh=calibration_tree.calibrationVectorList.calibrationVector(ii).line;
        if hh==hh1
            cali(b2,1:co)=calibration_tree.calibrationVectorList.calibrationVector(ii).sigmaNought.CONTENT;            
            hh1=calibration_tree.calibrationVectorList.calibrationVector(ii).line;           
        else
            b2=b2+1;           
            cali(b2,1:co)=calibration_tree.calibrationVectorList.calibrationVector(ii).sigmaNought.CONTENT;           
            hh1=calibration_tree.calibrationVectorList.calibrationVector(ii).line;
        end
    end
end
cali(sum(cali,1)==0,:)=[];
[rows,cols]=size(cali);
[cx,cy]=meshgrid(1:1:rows,1:1:cols);
rowss=xf1;colss=yf1;
[cxx,cyy]=meshgrid(1:((rows-1)/(rowss-1)):rows,1:((cols-1)/(colss-1)):cols);
cali1=interp2(cx,cy,cali',cxx,cyy,'cubic');
cali=cali1';
clear  b1 b2 ll ii hh hh1 co cxx cyy cx cy rows cols rowss colss cali1

%% read and interpolate thermal noise
b2=1; double no[]; 
no2=[];
if time < 20180313    
% If the IPF version priori to the version 2.9,only the noise vector in range 
% direction is provided by ESA
ll=length(noise_tree.noiseVectorList.noiseVector);
    for ii=1:ll    
        if ii==1
            co=noise_tree.noiseVectorList.noiseVector(ii).noiseLut.ATTRIBUTE.count;
            noo=noise_tree.noiseVectorList.noiseVector(ii).noiseLut.CONTENT;
            no(b2,1:co)=noo;
            no2(b2,1:co)=prctile(noo(noo(:)>0),5);
            hh1=noise_tree.noiseVectorList.noiseVector(ii).line;        
        else
            co=noise_tree.noiseVectorList.noiseVector(ii).noiseLut.ATTRIBUTE.count;
            hh=noise_tree.noiseVectorList.noiseVector(ii).line;
            noo=noise_tree.noiseVectorList.noiseVector(ii).noiseLut.CONTENT;  
            if hh==hh1
                no(b2,1:co)=noo;
                no2(b2,1:co)=prctile(noo(noo(:)>0),5);
                hh1=noise_tree.noiseVectorList.noiseVector(ii).line;            
            else
                b2=b2+1;           
               no(b2,1:co)=noo;       
               no2(b2,1:co)=prctile(noo(noo(:)>0),5);
               hh1=noise_tree.noiseVectorList.noiseVector(ii).line;           
            end
        end
    end
else
%For the IPF version after the version 2.9(available from Mar.13, 2018)
%The noise vectors in both range and azimuth direction are available.
    leng=length(noise_tree.noiseRangeVectorList.noiseRangeVector);
    RangeXn=zeros();RangeYn=zeros();sigma0n1=zeros();
    for ii=1:leng
        RangeXn(ii)=noise_tree.noiseRangeVectorList.noiseRangeVector(ii).line+1;
        RangeYnn=noise_tree.noiseRangeVectorList.noiseRangeVector(ii).pixel;
        RangeYn(ii,1:RangeYnn.ATTRIBUTE.count)=RangeYnn.CONTENT+1;
        sigma011=noise_tree.noiseRangeVectorList.noiseRangeVector(ii).noiseRangeLut;   
        sigma0n1(ii,1:RangeYnn.ATTRIBUTE.count)=sigma011.CONTENT;
    end 
    
    %process range noise vectors at the frames around the image
    [~,y]=find(sigma0n1==0);
    if isempty(y)==0                                %if noise vector has 0 value
        for i=2:length(y) 
            if y(i)~=y(i-1) 
                sigma0n1oneline=sigma0n1(:,y(i-1));
                index=find(sigma0n1oneline==0);     %find the position of 0 value and remove them
                sigma0n1oneline(index)=[];
                xx=RangeXn;
                xx2=xx(index);
                xx(index)=[];
                if length(sigma0n1oneline)>10
                    sigma0n1(index,y(i-1))=interp1(xx,sigma0n1oneline,xx2,'pchip');
                end
            end
        end
    end
   
    %interpolate range noise and azimuth noise
    clear('azNoise','azNoise11');
    leng=length(noise_tree.noiseAzimuthVectorList.noiseAzimuthVector);
    no=zeros();  no2=zeros();
    %calculate the complete noise in each sub-swath
    for ii=1:leng 
        clear('firstAzimuthLine','firstRangeSample','lastAzimuthLine','lastRangeSample');
        firstAzimuthLine=noise_tree.noiseAzimuthVectorList.noiseAzimuthVector(ii).firstAzimuthLine+1;
        firstRangeSample=noise_tree.noiseAzimuthVectorList.noiseAzimuthVector(ii).firstRangeSample+1;
        lastAzimuthLine=noise_tree.noiseAzimuthVectorList.noiseAzimuthVector(ii).lastAzimuthLine+1;
        lastRangeSample=noise_tree.noiseAzimuthVectorList.noiseAzimuthVector(ii).lastRangeSample+1; 
        clear('x','x1','x2','index1','index2','index11','index22','index3','index4');
        
        %find the closest azimuth line in calibrating matrix to the first azimuth line in subswath
        [x1,index11]=sort(abs(RangeXn-firstAzimuthLine));
        if RangeXn(index11(1))<=firstAzimuthLine
            index1=index11(1);
        else
            index1=index11(1)-1;
        end
        
        %find the closest azimuth line in calibrating matrix to the last azimuth line in subswath
        [x2,index22]=sort(abs(RangeXn-lastAzimuthLine));
        if RangeXn(index22(1))<lastAzimuthLine
            index2=index22(1)+1;
        else
            index2=index22(1);
        end
        
        %find the closest azimuth line in calibrating matrix to the last azimuth line in subswath
        if index2-index1<=1 
            if x1(1)<x2(1)
                midAzsubsw=index11(1);
            else
                midAzsubsw=index22(1);
            end
            clear('index3','index4','temp','rrsam');
            [~,index3]=sort(abs(RangeYn(midAzsubsw,:)-firstRangeSample));
            [~,index4]=sort(abs(RangeYn(midAzsubsw,:)-lastRangeSample));
            rrsam=firstRangeSample:1:lastRangeSample;
            temp=interp1(RangeYn(midAzsubsw,index3(1):index4(1)),sigma0n1(midAzsubsw,index3(1):index4(1)),rrsam,'pchip');
            for jj=1:lastAzimuthLine-firstAzimuthLine+1
               rsNoise1(jj,:)=temp; 
            end
        
        %if the azimuth width is equal to 1 pixel, then the range vector is the 1-D interpolation of the closest range vector
        elseif index2-index1==2 
            clear('index3','index4');
            midAzsubsw=ceil((index1+index2)/2);
            index3=find(RangeYn(midAzsubsw,:)==firstRangeSample);
            index4=find(RangeYn(midAzsubsw,:)==lastRangeSample);

            nx=RangeXn(index1:index2);ny=RangeYn(midAzsubsw,index3:index4);
            [nY,nX]=meshgrid(ny,nx);
            nxx=firstAzimuthLine:1:lastAzimuthLine;nyy=firstRangeSample:1:lastRangeSample;
            [nYY,nXX]=meshgrid(nyy,nxx);
            rsNoise1=griddata(nX,nY,sigma0n1(index1:index2,index3:index4),nXX,nYY,'cubic');
        
        %if the azimuth width is more than 1 pixel
        else
            clear('index3','index4','temp','sigma0n2','nx','ny','nxx','nyy','nX','nY','nXX','nYY');
            midAzsubsw=ceil((index1+index2)/2);
            index3=find(RangeYn(midAzsubsw,:)==firstRangeSample);
            index4=find(RangeYn(midAzsubsw,:)==lastRangeSample);
            sigma0n2=zeros(size(sigma0n1(index1:index2,index3:index4)));
            for mm=index3:index4
                nx=zeros(1,length(RangeXn(index1:index2)));
                nx(1,1)=firstAzimuthLine;
                nx(1,2:end-1)=RangeXn(index1+1:index2-1);
                nx(1,end)=lastAzimuthLine;
                temp=interp1(nx(1,2:end-1),sigma0n1(index1+1:index2-1,mm)',nx,'pchip');
                sigma0n2(:,mm-index3+1)=temp';
            end
            ny=RangeYn(midAzsubsw,index3:index4);
            [nY,nX]=meshgrid(ny,nx);
            nxx=firstAzimuthLine:1:lastAzimuthLine;nyy=firstRangeSample:1:lastRangeSample;
            [nYY,nXX]=meshgrid(nyy,nxx);
            rsNoise1=griddata(nX,nY,sigma0n2,nXX,nYY,'cubic');     
        end
        
        %interpolate the azimuth noise
        azLine=noise_tree.noiseAzimuthVectorList.noiseAzimuthVector(ii).line;
        noiseAzimuthLut=noise_tree.noiseAzimuthVectorList.noiseAzimuthVector(ii).noiseAzimuthLut;
        noiseAzimuthLutval=noiseAzimuthLut.CONTENT;
        if ischar(noiseAzimuthLutval)
            noiseAzimuthLutval2=zeros();
            index1=zeros(1,azLine.ATTRIBUTE.count+1);
            index1(2:end-1)=find(isspace(noiseAzimuthLutval));
            index1(1)=0;index1(end)=length(noiseAzimuthLutval)+1;
            for jj2=1:length(index1)-1
                noiseAzimuthLutval2(jj2)=str2double(noiseAzimuthLutval(index1(jj2)+1:index1(jj2+1)-1));
            end
            noiseAzimuthLutval=noiseAzimuthLutval2;
        end
        azline2=azLine.CONTENT;
        if ischar(azline2)
            azline2temp=zeros();
            index1=zeros(1,azLine.ATTRIBUTE.count+1);
            index1(2:end-1)=find(isspace(azline2));
            index1(1)=0;index1(end)=length(azline2)+1;
            for jj2=1:length(index1)-1
                azline2temp(jj2)=str2double(azline2(index1(jj2)+1:index1(jj2+1)-1));
            end
            azline2=azline2temp;
        end
        
        %if the rows number of subswath is equal to 1, then all the subswath value are same
        if azLine.ATTRIBUTE.count==1 
            azNoise1=noiseAzimuthLutval;
        else                                                              %the azimuth noise is 2-D interpolation of the azimuth noise vector
            azxx=firstAzimuthLine:1:lastAzimuthLine;
            azNoise11=interp1(azline2,noiseAzimuthLutval,azxx,'pchip');
            swathWide=lastRangeSample-firstRangeSample+1;
            azNoise22=repmat(azNoise11,swathWide,1);
            azNoise1=azNoise22';
        end
        
        noo=azNoise1.*rsNoise1;                                           %thermal noise is equal to azimuth noise multiply to range noise 
        no(firstAzimuthLine:lastAzimuthLine,firstRangeSample:lastRangeSample)=noo;
        no2(firstAzimuthLine:lastAzimuthLine,firstRangeSample:lastRangeSample)=prctile(noo(noo(:)>0),5);
        clear azNoise1 rsNoise1 noo
    end
    clear('RangeYnn','sigma011','midAzsubsw','rrsam','temp','nx','ny','nxx','nyy','nX','nY','nXX','nYY','azLine','noiseAzimuthLut','azNoise11','azNoise22','azxx','swathWide');
    clear('firstAzimuthLine','firstRangeSample','lastAzimuthLine','lastRangeSample');
    clear('x','x1','x2','index1','index2','index11','index22','index3','index4');
end
[rows,cols]=size(no);
if xf1~=rows||yf1~=cols
    [cx,cy]=meshgrid(1:rows,1:cols);
    rowss=xf1;colss=yf1;
    [cxx,cyy]=meshgrid(1:((rows-1)/(rowss-1)):rows,1:((cols-1)/(colss-1)):cols);
    no1=interp2(cx,cy,no',cxx,cyy,'cubic'); 
    no=no1';
    no2=interp2(cx,cy,no2',cxx,cyy,'cubic'); 
    no2=no2';
end
clear  b1 b2 ll ii hh hh1 co cxx cyy cx cy rows cols rowss colss no1 xf1 yf1 
if isempty(no2)
    no2=no;
end
% thermal noise removal
sigma0L_HH=(double(Image.^2)-no)./cali.^2;
sigma0L_HH(sigma0L_HH<0)=0;
noise=no2./cali.^2;
clear  no cali calibration_tree noise_tree


%% black edge removal 
[xf1,yf1]=size(Image);  
 for ii=1:yf1
    id=find(Image(:,ii)<=20);
    if length(id)/xf1>0.15
        idy=ii;
    else
        
        idy=ii;
        break;
    end
end

for ii=1:xf1
    id=find(Image(ii,:)<=20);
    if length(id)/yf1>0.15
        idx=ii;
    else   
        idx=ii;
        break;
    end
end
for ii=xf1:-1:1
     id=find(Image(ii,:)<=20);
    if length(id)/yf1>0.15
        idxh=ii;
    else
         idxh=ii;
        break;  
    end
end
for ii=yf1:-1:1
     id=find(Image(:,ii)<=20);
    if length(id)/xf1>0.15
        idyh=ii;
    else
        idyh=ii;
        break;
    end
end        
sigma0L_hh=sigma0L_HH(idx:idxh,idy:idyh);
noise_nesz=noise(idx:idxh,idy:idyh);
clear id  ii xf1 yf1

%% interpolate latitude, longitude and incident angle matrix
clear('lat','lon','lat2','lon2','rowline','colpixel');
  n=1;m=1;
lat=zeros;lon=zeros;rowline=zeros;colpixel=zeros;
ll=length(meta.geolocationGrid.geolocationGridPointList.geolocationGridPoint);
%read the latitude, longitude and incident angle vector
for ii=1:ll
    lat(n,m)=meta.geolocationGrid.geolocationGridPointList.geolocationGridPoint(ii).latitude;
    lon(n,m)=meta.geolocationGrid.geolocationGridPointList.geolocationGridPoint(ii).longitude;
    inc(n,m)=meta.geolocationGrid.geolocationGridPointList.geolocationGridPoint(ii).incidenceAngle;
    rowline(n,m)=meta.geolocationGrid.geolocationGridPointList.geolocationGridPoint(ii).line+1;
    colpixel(n,m)=meta.geolocationGrid.geolocationGridPointList.geolocationGridPoint(ii).pixel+1;
    if ii<ll
        if rowline(n,m)~=meta.geolocationGrid.geolocationGridPointList.geolocationGridPoint(ii+1).line+1
            n=n+1;
            m=1;
        else 
            m=m+1;
        end
    end
end

[xf1,yf1]=size(sigma0L_HH);
nxx=1:1:xf1;nyy=1:1:yf1;
[nY,nX]=meshgrid(nyy,nxx);
lat2=griddata(rowline,colpixel,lat,nX,nY,'cubic');                        %interpolate the latitude matrix

% if the longitude crosses 180 degree,then convert it to 0-360 degrees
if max(max(lon))>0 && min(min(lon))<0 && abs(max(max(lon))-180)<5 && abs(min(min(lon))-(-180)<5)
    lon=rem(360+lon,360);
    lon2=griddata(rowline,colpixel,lon,nX,nY,'cubic');                    %interpolate the longitude matrix
else
    lon2=griddata(rowline,colpixel,lon,nX,nY,'cubic');
    lon2=rem(360+lon2,360);  
end

% if the longitude crosses 0 degree,then convert it to -180-180 degrees
if max(max(lon2))>350 && min(min(lon2))<10
    lon2(lon2>180)=lon2(lon2>180)-360;  
end

%interpolate incident angle matrix
[rows,cols]=size(inc);
[cx,cy]=meshgrid(1:1:rows,1:1:cols);
[cxx,cyy]=meshgrid(1:((rows-1)/(xf1-1)):rows,1:((cols-1)/(yf1-1)):cols);
inc1=interp2(cx,cy,inc',cxx,cyy,'cubic');
inc2=inc1';

%remove the black edge of latitude, longitude and incident angle
lat2=lat2(idx:idxh,idy:idyh);
lon2=lon2(idx:idxh,idy:idyh);
inc2=inc2(idx:idxh,idy:idyh);

clear rowline colpixel xf1 yf1 nX nY nxx nyy lon lat inc inc1 idx idxh idy idyh rows cols cx cxx cy cyy
end

