function ctSent_deburst
%CTSENT_DEBURST Summary of this function goes here
%   Detailed explanation goes here
fprintf('#######################CtSent v1.0####################### \n');
fprintf('######################################################### \n');
fprintf('###################     Deburst    ###################### \n');
fprintf('######################################################### \n');
fprintf('########     Zhang-Feng Ma, Hohai University      ####### \n');
fprintf('########        Mi Jiang  , Hohai University      ####### \n');
fprintf('########         Yi Zhao  , Hohai University      ####### \n');
fprintf('###################  30,July,2019    #################### \n');
parm = load('parms');
filepath = load('intfpath');
slcpath = filepath.intfpath;
geopath = load('geopath');
geopath = geopath.geopath;
intfpath = load('intfpath');
intfpath = intfpath.intfpath;
ndate = length(slcpath);
masterIndex = 0;
for j = 1:ndate
    datestr = slcpath(j).date;
    if strcmp(datestr,parm.masterdate)
        masterIndex = j;
    end
end
if masterIndex == 0 
    fprintf('Master date cannot be detected and set the first acquisition as the Master reference... \n');
    masterIndex = 1;
end
infopath = load('slcpath');
infopath = infopath.slcpath;
info_mas = load(infopath(masterIndex).infopath);
info_mas = info_mas.infosave;
[ Lines,Pixels,Rect] = TOPSmerge(info_mas);
nsubswath = length(info_mas.IW);
fileDir  = [parm.workDir,'\InSAR',info_mas.date,'\geo'];
if ~exist(fileDir,'dir')
   mkdir(fileDir);
end
debGeo.lat = [fileDir,'\lat'];
debGeo.lon = [fileDir,'\lon'];
debGeo.hgt = [fileDir,'\hgt'];
debGeo.inc = [fileDir,'\inc'];
LatMerge = zeros(Lines,Pixels,'single');
for j = 1:nsubswath
    lat = [];
    for k = 1:info_mas.IW(j).subSwath.numOfBursts
      data = freadCt(geopath.Lat{j,k},info_mas.IW(j).sensor_para.linesPerBurst,info_mas.IW(j).sensor_para.samplesPerBurst);
      lat = [lat;data];
    end
    latdeb = TOPSAR_DEBURST(info_mas.IW(j),1,info_mas.IW(j).subSwath.numOfBursts,lat);
    LatMerge(Rect(j).RectInOutBand(1):Rect(j).RectInOutBand(3),Rect(j).RectInOutBand(2):Rect(j).RectInOutBand(4)) = ...
        latdeb(Rect(j).RectInInput(1):Rect(j).RectInInput(3),Rect(j).RectInInput(2):Rect(j).RectInInput(4));
end
% fwriteCt(debGeo.lat,LatMerge);
clear lat;clear data;clear latdeb;
% clear LatMerge
fclose('all');
LonMerge = zeros(Lines,Pixels,'single');
for j = 1:nsubswath
    lon = [];
    for k = 1:info_mas.IW(j).subSwath.numOfBursts
      data = freadCt(geopath.Lon{j,k},info_mas.IW(j).sensor_para.linesPerBurst,info_mas.IW(j).sensor_para.samplesPerBurst);
      lon = [lon;data];
    end
    londeb = TOPSAR_DEBURST(info_mas.IW(j),1,info_mas.IW(j).subSwath.numOfBursts,lon);
    LonMerge(Rect(j).RectInOutBand(1):Rect(j).RectInOutBand(3),Rect(j).RectInOutBand(2):Rect(j).RectInOutBand(4)) = ...
        londeb(Rect(j).RectInInput(1):Rect(j).RectInInput(3),Rect(j).RectInInput(2):Rect(j).RectInInput(4));
end
% fwriteCt(debGeo.lon,LonMerge);
clear lon;clear data;clear londeb;
% clear LonMerge
fclose('all');
[ROI] = getExtractROI(LatMerge,LonMerge,parm.rectangle);
fwriteCt(debGeo.lat,LatMerge(ROI(1):ROI(2),ROI(3):ROI(4)));
fwriteCt(debGeo.lon,LonMerge(ROI(1):ROI(2),ROI(3):ROI(4)));
clear LatMerge;
clear LonMerge;
HgtMerge = zeros(Lines,Pixels,'single');
for j = 1:nsubswath
    hgt = [];
    for k = 1:info_mas.IW(j).subSwath.numOfBursts
      data = freadCt(geopath.Hgt{j,k},info_mas.IW(j).sensor_para.linesPerBurst,info_mas.IW(j).sensor_para.samplesPerBurst);
      hgt= [hgt;data];
    end
    hgtdeb = TOPSAR_DEBURST(info_mas.IW(j),1,info_mas.IW(j).subSwath.numOfBursts,hgt);
    HgtMerge(Rect(j).RectInOutBand(1):Rect(j).RectInOutBand(3),Rect(j).RectInOutBand(2):Rect(j).RectInOutBand(4)) = ...
        hgtdeb(Rect(j).RectInInput(1):Rect(j).RectInInput(3),Rect(j).RectInInput(2):Rect(j).RectInInput(4));
end
fwriteCt(debGeo.hgt,HgtMerge(ROI(1):ROI(2),ROI(3):ROI(4)));
clear hgt;clear data;clear hgtdeb;clear HgtMerge
IncMerge = zeros(Lines,Pixels,'single');
for j = 1:nsubswath
    inc = [];
    for k = 1:info_mas.IW(j).subSwath.numOfBursts
      data = freadCt(geopath.Inc{j,k},info_mas.IW(j).sensor_para.linesPerBurst,info_mas.IW(j).sensor_para.samplesPerBurst);
      inc= [inc;data];
    end
    incdeb = TOPSAR_DEBURST(info_mas.IW(j),1,info_mas.IW(j).subSwath.numOfBursts,inc);
    IncMerge(Rect(j).RectInOutBand(1):Rect(j).RectInOutBand(3),Rect(j).RectInOutBand(2):Rect(j).RectInOutBand(4)) = ...
        incdeb(Rect(j).RectInInput(1):Rect(j).RectInInput(3),Rect(j).RectInInput(2):Rect(j).RectInInput(4));
end
fwriteCt(debGeo.inc,IncMerge(ROI(1):ROI(2),ROI(3):ROI(4)));
clear inc;clear data;clear incdeb;clear IncMerge
fclose('all');
fileDir  = [parm.workDir,'\InSAR',info_mas.date,'\slc'];
if ~exist(fileDir,'dir')
   mkdir(fileDir);
end
for i = 1:ndate
SlcMerge = zeros(Lines,Pixels,'single');
slcsavepath(i).savepath = [fileDir,'\',slcpath(i).date,'.slc'];
for j = 1:nsubswath
    slc = [];
    for k = 1:info_mas.IW(j).subSwath.numOfBursts
        slaveBurstIndex = info_mas.IW(j).BurstUse(i,k);
        if slaveBurstIndex
        data = readcpx(slcpath(i).savepath{j,k},info_mas.IW(j).sensor_para.samplesPerBurst);
        else
        data = zeros(info_mas.IW(j).sensor_para.linesPerBurst,info_mas.IW(j).sensor_para.samplesPerBurst,'single');
        end
      slc= [slc;data];
    end
    slcdeb = TOPSAR_DEBURST(info_mas.IW(j),1,info_mas.IW(j).subSwath.numOfBursts,slc);
    clear slc;
    SlcMerge(Rect(j).RectInOutBand(1):Rect(j).RectInOutBand(3),Rect(j).RectInOutBand(2):Rect(j).RectInOutBand(4)) = ...
        slcdeb(Rect(j).RectInInput(1):Rect(j).RectInInput(3),Rect(j).RectInInput(2):Rect(j).RectInInput(4));
end
writecpx(slcsavepath(i).savepath,SlcMerge(ROI(1):ROI(2),ROI(3):ROI(4)));
clear slcdeb;clear SlcMerge
fclose('all');
end
fileDir  = [parm.workDir,'\InSAR',info_mas.date,'\baseline'];
if ~exist(fileDir,'dir')
   mkdir(fileDir);
end
for i = 1:ndate
BaseMerge = zeros(Lines,Pixels,'single');
basesavepath(i).savepath = [fileDir,'\',slcpath(i).date,'.base'];
for j = 1:nsubswath
    base = [];
    for k = 1:info_mas.IW(j).subSwath.numOfBursts
    slaveBurstIndex = info_mas.IW(j).BurstUse(i,k);
      if slaveBurstIndex
        data = freadCt(intfpath(i).basepath{j,k},info_mas.IW(j).sensor_para.linesPerBurst,info_mas.IW(j).sensor_para.samplesPerBurst);
      else
        data = zeros(info_mas.IW(j).sensor_para.linesPerBurst,info_mas.IW(j).sensor_para.samplesPerBurst,'single');
      end
      base = [base;data];
    end
    basedeb = TOPSAR_DEBURST(info_mas.IW(j),1,info_mas.IW(j).subSwath.numOfBursts,base);
    clear base;
    BaseMerge(Rect(j).RectInOutBand(1):Rect(j).RectInOutBand(3),Rect(j).RectInOutBand(2):Rect(j).RectInOutBand(4)) = ...
        basedeb(Rect(j).RectInInput(1):Rect(j).RectInInput(3),Rect(j).RectInInput(2):Rect(j).RectInInput(4));
end
fwriteCt(basesavepath(i).savepath,BaseMerge(ROI(1):ROI(2),ROI(3):ROI(4)));
clear basedeb;clear BaseMerge
fclose('all');
end
Lines = ROI(2)-ROI(1)+1;
Pixels = ROI(4)-ROI(3)+1;
save 'debGeo.mat' debGeo Lines Pixels;
save 'debSlc.mat' slcsavepath Lines Pixels;
save 'baseline.mat' basesavepath Lines Pixels;
end
function [imageDeb] = TOPSAR_DEBURST( info,BurstStart,BurstEnd,input)
tic;
starty0=(BurstStart-1)*info.subSwath.linesPerBurst;
y0=round((info.subSwath.burstFirstLineTime(BurstStart)-info.subSwath.burstFirstLineTime(1))/info.subSwath.azimuthTimeInterval)+1;
yMax=round((info.subSwath.burstLastLineTime(BurstEnd)-info.subSwath.burstFirstLineTime(1))/info.subSwath.azimuthTimeInterval)+1;
nlines = yMax - y0 + 1;
[~,~,npages] = size(input);
imageDeb = zeros(nlines,info.sensor_para.samplesPerBurst,npages,'single');
for j=y0:yMax
 targetLineTime=info.subSwath.burstFirstLineTime(1)+(j-1)*info.subSwath.azimuthTimeInterval;
 burstTimes.sy0 = -1;
 burstTimes.sy1 = -1;
k = 0;
        for i = BurstStart:BurstEnd
            if (targetLineTime >= (info.subSwath.burstFirstLineTime(i)-0.00021) && targetLineTime < (info.subSwath.burstLastLineTime(i)+0.00021))
                sy = (i-1) * info.subSwath.linesPerBurst +...
                        floor(((targetLineTime - info.subSwath.burstFirstLineTime(i)) / info.subSwath.azimuthTimeInterval)+0.5);
                if (k == 0) 
                    burstTimes.sy0 = sy;
                    burstTimes.burstNum0 = i;
                 else 
                    burstTimes.sy1 = sy;
                    burstTimes.burstNum1 = i;
                    break;
                end
                k=k+1;
            end
        end
        if (burstTimes.sy0 == -1 && burstTimes.sy1 == -1) 
            continue;
        end
        if (burstTimes.sy0 ~= -1 && burstTimes.sy1 ~= -1) 
            burstTimes.midTime = (info.subSwath.burstLastLineTime(burstTimes.burstNum0) +...
                    info.subSwath.burstFirstLineTime(burstTimes.burstNum1)) / 2.0;
        end
        if (burstTimes.sy1 ~= -1 && targetLineTime > burstTimes.midTime) 
                imageDeb(j-y0+1,:) = input(burstTimes.sy1+1-starty0,:,:);
             else 
                imageDeb(j-y0+1,:) = input(burstTimes.sy0+1-starty0,:,:);
        end

end
t=toc;
        fprintf('TOPS Deburst: %d seconds finished...\n',floor(t));
end
function fwriteCt(savepath,savedata)
fileID = fopen(savepath,'w');
fwrite(fileID,savedata,'float');
fclose(fileID);
end
function A = freadCt(savepath,lines,width)
fileID = fopen(savepath);
A = fread(fileID,[lines width],['float','=>single']);
fclose(fileID);
end
function [ Lines,Pixels,Rect] = TOPSmerge(info)
%Zhangfeng Ma,Mi Jiang 2018/9/7
%Hohai University, School of Earth Science and Engineering
%E-Mail;jspcmazhangfeng@hhu.edu.cn;mijiang@hhu.edu.cn
[Lines,Pixels,RangeStartTime,RangeEndTime,first_line_time,...
    last_line_time,first_pixel_time,last_pixel_time] = getOutBandSize(info);

[Rect] = getCorrespondREC(info,RangeStartTime,RangeEndTime,first_line_time,...
    first_pixel_time);
end
function [Lines,Pixels,RangeStartTime,RangeEndTime,first_line_time,...
    last_line_time,first_pixel_time,last_pixel_time] = getOutBandSize(info)
nsubswath = length(info.IW);
FirstRangeTime = zeros(nsubswath,1);
LastRangeTime = zeros(nsubswath,1);
FirstLineTime = zeros(nsubswath,1);
LastLineTime = zeros(nsubswath,1);
for i = 1:nsubswath
    FirstRangeTime(i) = info.IW(i).subSwath.slrTimeToFirstPixel;
    LastRangeTime(i) = info.IW(i).subSwath.slrTimeToLastPixel;
    FirstLineTime(i)= info.IW(i).subSwath.burstFirstLineTime(1);
    LastLineTime(i) = info.IW(i).subSwath.burstLastLineTime(end);
end
first_pixel_time = min(FirstRangeTime);
last_pixel_time = max(LastRangeTime);
rangetimespacing = info.IW(1).subSwath.rangePixelSpacing/299792458;
Pixels = round((last_pixel_time - first_pixel_time)/rangetimespacing) + 1;
first_line_time = min(FirstLineTime);
last_line_time = max(LastLineTime);
Lines = round((last_line_time - first_line_time)/...
    info.IW(1).sensor_para.line_time_interval) + 1;
RangeStartTime = zeros(nsubswath,1);
RangeEndTime = zeros(nsubswath,1);
if nsubswath > 1
for i = 1:nsubswath-1
    if i == 1
        RangeStartTime(i) = FirstRangeTime(i);
    end
    if i == nsubswath-1
        RangeEndTime(i+1) = LastRangeTime(i+1);
    end
   LapNum = round((LastRangeTime(i) - FirstRangeTime(i+1))/rangetimespacing);
   halfLapNum = floor(LapNum/2);
   RangeEndTime(i) = LastRangeTime(i) - halfLapNum * rangetimespacing;
   RangeStartTime(i+1) = FirstRangeTime(i+1) + halfLapNum * rangetimespacing;
end
else
    RangeEndTime = LastRangeTime;
    RangeStartTime = FirstRangeTime;
end
end
function [Rect] = getCorrespondREC(info,RangeStartTime,RangeEndTime,first_line_time,...
    first_pixel_time)
for i = 1:length(info.IW)
rangetimespacing = info.IW(1).subSwath.rangePixelSpacing/299792458;
tAziStart = info.IW(i).subSwath.burstFirstLineTime(1);
tAziEnd = info.IW(i).subSwath.burstLastLineTime(end);
 
%%%%%%%%%%%%%%%%  Rectangle in OUTBAND  %%%%%%%%%%%%%%%%%%%
yOUT0 = round((tAziStart-first_line_time)/info.IW(1).sensor_para.line_time_interval+1);
xOUT0 = round((RangeStartTime(i)-first_pixel_time)/rangetimespacing+1);
yOUT1 = round((tAziEnd-first_line_time)/info.IW(1).sensor_para.line_time_interval+1);
xOUT1 = round((RangeEndTime(i)-first_pixel_time)/rangetimespacing+1);
Rect(i).RectInOutBand = [yOUT0,xOUT0,yOUT1,xOUT1];
%%%%%%%%%%%%%%%% Rectangle in INPUTIMAGE %%%%%%%%%%%%%%%%%%%%
yIN0 = 1;
xIN0 = round((RangeStartTime(i)-info.IW(i).subSwath.slrTimeToFirstPixel)/rangetimespacing+1);
yIN1 = round((tAziEnd-tAziStart)/info.IW(1).sensor_para.line_time_interval+1);
xIN1 = round((RangeEndTime(i)-info.IW(i).subSwath.slrTimeToFirstPixel)/rangetimespacing+1);
Rect(i).RectInInput = [yIN0,xIN0,yIN1,xIN1];
end
end
function [ROI] = getExtractROI(latmerge,lonmerge,rectangle)
lat0 = rectangle(1);lat1 = rectangle(2);
lon0 = rectangle(3);lon1 = rectangle(4);
[nline,npixel] = size(latmerge);
azimuthline = (1:5:nline)';
rangepixel = 1:20:npixel;
Line = repmat(azimuthline,1,length(rangepixel));
Pixel = repmat(rangepixel,length(azimuthline),1);
leftlowerPixelID = round(griddata(double(latmerge(1:5:end,1:20:end)),double(lonmerge(1:5:end,1:20:end)),Pixel,lat0,lon0));
leftlowerLineID = round(griddata(double(latmerge(1:5:end,1:20:end)),double(lonmerge(1:5:end,1:20:end)),Line,lat0,lon0));
leftupperPixelID = round(griddata(double(latmerge(1:5:end,1:20:end)),double(lonmerge(1:5:end,1:20:end)),Pixel,lat1,lon0));
leftupperLineID = round(griddata(double(latmerge(1:5:end,1:20:end)),double(lonmerge(1:5:end,1:20:end)),Line,lat1,lon0));
rightlowerPixelID = round(griddata(double(latmerge(1:5:end,1:20:end)),double(lonmerge(1:5:end,1:20:end)),Pixel,lat0,lon1));
rightlowerLineID = round(griddata(double(latmerge(1:5:end,1:20:end)),double(lonmerge(1:5:end,1:20:end)),Line,lat0,lon1));
rightupperPixelID = round(griddata(double(latmerge(1:5:end,1:20:end)),double(lonmerge(1:5:end,1:20:end)),Pixel,lat1,lon1));
rightupperLineID = round(griddata(double(latmerge(1:5:end,1:20:end)),double(lonmerge(1:5:end,1:20:end)),Line,lat1,lon1));
if leftupperLineID > leftlowerLineID
    lineID1 = leftupperLineID;
    lineID0 = rightlowerLineID;
    pixelID1 = rightupperPixelID;
    pixelID0 = leftlowerPixelID;
else
    lineID1 = leftlowerLineID;
    lineID0 = rightupperLineID;
    pixelID1 = leftupperPixelID;
    pixelID0 = rightlowerPixelID;
end
ROI = [lineID0,lineID1,pixelID0,pixelID1];
end
function [vq] = trinterp(x,y,v,xq,yq,ratio)
[m,n] = size(xq);
X = x(:);Y = y(:)*ratio;Z = v(:);
xq = xq(:);yq = yq(:);
TRI = delaunayTriangulation(X,Y);
ID = pointLocation(TRI,xq,yq.*ratio);
Vx0 = X(TRI(ID,1));Vx1 = X(TRI(ID,2));Vx2 = X(TRI(ID,3));
Vy0 = Y(TRI(ID,1))./ratio;Vy1 = Y(TRI(ID,2))./ratio;Vy2 = Y(TRI(ID,3))./ratio;
Vz0 = Z(TRI(ID,1));Vz1 = Z(TRI(ID,2));Vz2 = Z(TRI(ID,3));
xkj = Vx1 - Vx0;
ykj = Vy1 - Vy0;
xlj = Vx2 - Vx0;
ylj = Vy2 - Vy0;
f = 1 ./ (xkj.*ylj - ykj.*xlj);
zj = Vz0;zk = Vz1;zl = Vz2;
zkj = zk - zj;
zlj = zl - zj;
a = -f.*(ykj.*zlj - zkj.*ylj);
b = -f.*(zkj.*xlj - xkj.*zlj);
c = -a.*Vx1 - b.*Vy1 + zk;
zout = a.* xq + b.* yq + c;
vq = reshape(zout,[m,n]);
end