function  offset = ctSent_EnhancedSpectralDiversity
fprintf('#######################CtSent v1.0####################### \n');
fprintf('######################################################### \n');
fprintf('###########     Enhanced Spectral Diversity      ######## \n');
fprintf('######################################################### \n');
fprintf('########     Zhang-Feng Ma, Hohai University      ####### \n');
fprintf('########        Mi Jiang  , Hohai University      ####### \n');
fprintf('########         Yi Zhao  , Hohai University      ####### \n');
fprintf('###################   23,July,2019    ################### \n');
parm = load('parms');
filepath = load('coregpath');
slcpath = filepath.slcpath;
ndate = length(slcpath);
masterIndex = 0;
for i = 1:ndate
    datestr = slcpath(i).date;
    if strcmp(datestr,parm.masterdate)
        masterIndex = i;
    end
end
fprintf('    Extracting overlapping region... \n');
info_mas = load(slcpath(masterIndex).infopath);
info_mas = info_mas.infosave;
nSubswath = length(info_mas.IW);
for i = 1:nSubswath
    BurstUse = info_mas.IW(i).BurstUse;
    nBurst = info_mas.IW(i).subSwath.numOfBursts;
    prodUse = prod(BurstUse);
    ovl= false(nBurst-1,1);
    idx = 1;
    for j = 1:nBurst-1
       ovl(j) = logical(prodUse(j)*prodUse(j+1));
        if sum(ovl)==0
        fprintf('Some Image cracked...\n');
        elseif ovl(j) ~= 0
        info_mas.IW(i).OvlUse(idx) = j;
        idx = idx + 1;
        end
    end  
end
for k = 1:nSubswath
for i = 1
%     :length(info_mas.IW(k).OvlUse)
    if info_mas.IW(k).OvlUse(i)
    [Rect_up,Rect_lo] = getOverlappedRectangles(info_mas.IW(k).OvlUse(i),info_mas.IW(k).subSwath);
    lines = info_mas.IW(k).subSwath.linesPerBurst;
    width = info_mas.IW(k).subSwath.samplesPerBurst;
    fprintf(['    Extracting ',info_mas.IW(k).Swathname,' No. %d overlapping region... \n'],i);
    for j = 1:ndate
        info = load(slcpath(j).infopath);
        info = info.infosave;
        ESD_Dir= [parm.workDir,'\ESD\',info.date,'\',info_mas.IW(k).Swathname];
        if ~exist(ESD_Dir,'dir')
        mkdir(ESD_Dir);
        end
        fprintf(['------------' info.date '------------ \n']);
        slcpath(j).ovlpath{k,i,1} = [ESD_Dir,'\Overlapping',num2str(info_mas.IW(k).OvlUse(i)),'.up'];
        slcpath(j).ovlpath{k,i,2} = [ESD_Dir,'\Overlapping',num2str(info_mas.IW(k).OvlUse(i)),'.lo'];
%         if exist(slcpath(j).ovlpath{k,i,1},'file')&exist(slcpath(j).ovlpath{k,i,2},'file')
%             continue;
%         end
        [ovl_up, ~] = freadbk(slcpath(j).savepath{k,info_mas.IW(k).OvlUse(i)},lines,'mph',...
            Rect_up(2),Rect_up(2)+Rect_up(4)-1,Rect_up(1),Rect_up(1)+Rect_up(3)-1);
        [ovl_lo, ~] = freadbk(slcpath(j).savepath{k,info_mas.IW(k).OvlUse(i)+1},lines,'mph',...
            Rect_lo(2),Rect_lo(2)+Rect_lo(4)-1,Rect_lo(1),Rect_lo(1)+Rect_lo(3)-1);
%         slcpath(j).ovlpath{k,i,1} = [ESD_Dir,'\Overlapping',num2str(info_mas.IW(k).OvlUse(i)),'.up'];
%         slcpath(j).ovlpath{k,i,2} = [ESD_Dir,'\Overlapping',num2str(info_mas.IW(k).OvlUse(i)),'.lo'];
        writecpx(slcpath(j).ovlpath{k,i,1},ovl_up);
        writecpx(slcpath(j).ovlpath{k,i,2},ovl_lo);
        fprintf('    Writing finished ... \n');
        fclose('all');
    end
    end
end
end
save 'ovlpath.mat' slcpath;
%########### Here is the network needed to be replaced###########%
[order] = sequential_netwrok(ndate,parm.NetworkNum);
[offset] = estimateAzoff_network(info_mas,slcpath,order,masterIndex);
%########### Here is the network needed to be replaced###########%
% save 'Azioffset.mat' offset;
end
function [offset] = estimateAzoff_network(info_mas,slcpath,order,masterIndex)
parm = load('parms');
ndate = length(slcpath);
azioff = zeros(ndate,length(info_mas.IW));
Pweight = zeros(ndate,length(info_mas.IW));
for i = 1:length(info_mas.IW)
%     novl = length(info_mas.IW(i).OvlUse);
     novl = 1;
    [n,~]= size(order);
    G_mat = zeros(n,ndate); d_mat = zeros(n,novl);
    for j = 1:novl
        for k = 1:n
    [Rect_up,Rect_lo] = getOverlappedRectangles(info_mas.IW(i).OvlUse(j),info_mas.IW(i).subSwath);
    a=order(k,1);b = order(k,2);
    info_m = load(slcpath(a).infopath);
    info_m = info_m.infosave.IW(i);
    spectralSeparation = computeSpectralSeparation(info_m);
    mburst_ovl_up = readcpx(slcpath(a).ovlpath{i,j,1},Rect_up(3),Rect_up(4)); 
    mburst_ovl_lo = readcpx(slcpath(a).ovlpath{i,j,2},Rect_lo(3),Rect_lo(4)); 
    sburst_ovl_up = readcpx(slcpath(b).ovlpath{i,j,1},Rect_up(3),Rect_up(4)); 
    sburst_ovl_lo = readcpx(slcpath(b).ovlpath{i,j,2},Rect_lo(3),Rect_lo(4));
    fclose('all');
    forDiff = earlymultilooking(mburst_ovl_up,sburst_ovl_up); 
    backDiff = earlymultilooking(mburst_ovl_lo,sburst_ovl_lo);
    doubleDifference = forDiff.*conj(backDiff);
    %Coherence Estimator needed to be replaced
    Coh = BoxcarCoherence(mburst_ovl_up,sburst_ovl_up);
    setlle= angle(doubleDifference);
    residual = setlle(setlle ~= 0 & Coh > parm.esd_coh_threshold);
    misrad = fitMisregistration(residual);
%     [misrad,pcoh] = periodogram(residual);
    Pweight(k,j) = 1;
    fprintf('    %d Periodogram Completed ..,\n',k);
    d_mat(k,j) = misrad /(2*pi*spectralSeparation*info_m.sensor_para.line_time_interval);
    G_mat(k,a) = -1;G_mat(k,b) = 1;
        end
    end
   Indice = 1:ndate;
   Indice(masterIndex) = [];
   d_mat = mean(d_mat,2);
   Pweight = mean(Pweight,2);
   Pweight = diag(Pweight);
   G_mat(:,masterIndex) = [];
%    off = IRLS(G_mat,d_mat);
   off = pinv(G_mat'*Pweight*G_mat)*G_mat'*Pweight*d_mat;
   offset = zeros(ndate,1);
   offset(Indice) = off;
   azioff(:,i) = offset;
end
correctAzshift(-azioff,slcpath,masterIndex);
end
function correctAzshift(offset,slcpath,masterIndex)

ndate = length(slcpath);
info_mas = load(slcpath(masterIndex).infopath);
info_mas = info_mas.infosave;
nsubswath = length(info_mas.IW);
for i = 1:ndate
   for j = 1:nsubswath
        nburst = info_mas.IW(j).subSwath.numOfBursts;
        novl = nburst - 1;
     if novl == 0
       fprintf('Only one burst, Enhanced Spectral Diversity skipped...\n');
     end
    for k = 1:nburst
        if i == masterIndex
            continue;
        end
        slaveBurstIndex = info_mas.IW(j).BurstUse(i,k);
        if slaveBurstIndex
        fprintf(['    Correcting Azimuth shift of ',info_mas.IW(j).Swathname,' Burst %d in No. %d image... \n'],k,i);
        lines = info_mas.IW(j).subSwath.linesPerBurst;
        width = info_mas.IW(j).subSwath.samplesPerBurst;
%         temp = strcat(slcpath(i).savepath{j,k},'.esd');
%         if exist(temp,'file')
%             continue;
%         end
%         slavecplx = readcpx(slcpath(i).savepath{j,k},width,lines);
        [slavecplx,~] = freadbk(slcpath(i).savepath{j,k},lines,'mph');
        fileID = fopen(slcpath(i).derampdemodpath{j,k});
        derampDemodPhase = fread(fileID,[lines,width],['float','=>single']);
        fclose(fileID);
        slcpath(i).savepath{j,k} = strcat(slcpath(i).savepath{j,k},'.esd');
%         slcpath(i).savepath{j,k} = slcpath(i).savepath{j,k};
        [ correctIMG ] = performOvlAzshift(offset(i,j),slavecplx,real(derampDemodPhase) );
        writecpx(slcpath(i).savepath{j,k},correctIMG);
        fprintf('    Corrected... \n');
        end
    end
    end
end
        save 'esdpath.mat' slcpath;
end
function [ correctIMG ] = performOvlAzshift(azOffset,slc,derampDemodPhase )
derampDemod = exp(1j*derampDemodPhase);
[height,width] = size(derampDemod);
derampDemodcplx = slc .*derampDemod;
phaseMat = computeShiftPhaseArray(azOffset,height,width);
col2 = complex(derampDemodPhase,zeros(size(derampDemod),'single'));
col1 = ifft(fft(derampDemodcplx) .* phaseMat);
col2 = ifft(fft(col2) .* phaseMat);
phase = exp(1j*(real(col2)));
correctIMG = col1 .* conj(phase);
% correctIMG = col1;
end
function phaseMat = computeShiftPhaseArray(shift,height,width)
phase = -2.0 * pi * shift /height;
halfHeight = floor(height * 0.5 + 0.5);
phaseK = zeros(height,1,'single');
        for k = 0:height-1
            if (k < halfHeight) 
                phaseK(k+1) = phase * k;
             else 
                phaseK(k+1) = phase * (k - height);
            end
        end
        phaseArray = complex(cos(phaseK),sin(phaseK));
        phaseMat = repmat(phaseArray,1,width);
end
function [x] = IRLS(G,b)
maxiter = 20;
[n,~] = size(b);
temp = ones(n,1);
p=diag(temp);
x=inv(G'*p*G)*G'*p*b;
for iter = 1:maxiter
%x=inv(G'*p*G)*G'*p*b;
temp1 = abs(G*inv(G'*p*G)*G'*p*b-b);
temp1(find(temp1==0))=1;
betta = 1./temp1;
p = diag(betta);
solution = inv(G'*p*G)*G'*p*b;
if max(abs(solution-x)) < 1e-8
    break;
end
x = solution;
end
end
function  misrad = fitMisregistration(x)
[counts,centers] = hist(double(x), 120); 
[mu,~] = fit(centers',counts','Gauss2');
x1 = -pi:0.01:pi; 
a1 = mu.a1;b1 = mu.b1;c1 = mu.c1;a2 = mu.a2;b2 = mu.b2;c2 = mu.c2;
y1 = a1 * exp(-((x1-b1)/c1).^2) + a2 * exp(-((x1-b2)/c2).^2);
misrad = x1(find(y1 == max(y1)));
end
function Coh = BoxcarCoherence(mas,sla)
    h = fspecial('gaussian',[11 11],1); 
    % gaussian weight with sigma=2
    % intensity image of m and s
    nu = filter2(h,mas.*conj(sla));
        % denominator
    de1 = filter2(h,mas.*conj(mas));
    de2 = filter2(h,sla.*conj(sla));
    cohcplx = nu./sqrt(de1.*de2);
    Coh = abs(cohcplx);
end
function [Rectangle1,Rectangle2]=getOverlappedRectangles(overlapIndex,subSwath)
firstValidPixelOfBurstOne = getBurstFirstValidPixel(overlapIndex,subSwath);
lastValidPixelOfBurstOne = getBurstLastValidPixel(overlapIndex,subSwath);
firstValidPixelOfBurstTwo = getBurstFirstValidPixel(overlapIndex + 1,subSwath);
lastValidPixelOfBurstTwo = getBurstLastValidPixel(overlapIndex + 1,subSwath);
firstValidPixel = max(firstValidPixelOfBurstOne, firstValidPixelOfBurstTwo);
lastValidPixel = min(lastValidPixelOfBurstOne, lastValidPixelOfBurstTwo);
x0 = firstValidPixel;
w = lastValidPixel - firstValidPixel + 1;
numOfInvalidLinesInBurstOne = subSwath.linesPerBurst -...
                subSwath.lastValidLine(overlapIndex) - 1;
numOfInvalidLinesInBurstTwo = subSwath.firstValidLine(overlapIndex + 1);
numOverlappedLines = computeBurstOverlapSize(overlapIndex,subSwath);
h = numOverlappedLines - numOfInvalidLinesInBurstOne - numOfInvalidLinesInBurstTwo;
% y0BurstOne =subSwath.linesPerBurst * (overlapIndex) - numOfInvalidLinesInBurstOne - h;
y0BurstOne =subSwath.linesPerBurst - numOfInvalidLinesInBurstOne - h;
y0BurstTwo =  numOfInvalidLinesInBurstTwo;
Rectangle1=[x0, y0BurstOne, w, h];
Rectangle2=[x0, y0BurstTwo, w, h];
end
function index=getBurstFirstValidPixel(burstIndex,subSwath)
[~,n]=size(subSwath.firstValidSample);
        for lineIdx = 1:n
            if (subSwath.firstValidSample(burstIndex,lineIdx) ~= -1) 
                index=subSwath.firstValidSample(burstIndex,lineIdx);
                break;
            end
        end
end
function index=getBurstLastValidPixel(burstIndex,subSwath)
[~,n]=size(subSwath.lastValidSample);
        for lineIdx = 1:n
            if (subSwath.lastValidSample(burstIndex,lineIdx) ~= -1) 
                index=subSwath.lastValidSample(burstIndex,lineIdx);
                break;
            end
        end
end
function numOverlappedLines = computeBurstOverlapSize(overlapIndex,subSwath)
endTime = subSwath.burstLastLineTime(overlapIndex);
startTime = subSwath.burstFirstLineTime(overlapIndex + 1);
numOverlappedLines =floor((endTime - startTime) / subSwath.azimuthTimeInterval);
end
function SpectralSeparation=computeSpectralSeparation(info)
tCycle = info.subSwath.linesPerBurst * info.subSwath.azimuthTimeInterval;
sumSpectralSeparation = 0.0;
for b = 1:(info.subSwath.numOfBursts)
    dopplerRate = computeDopplerRate( info,b);
    sumSpectralSeparation = sumSpectralSeparation + dopplerRate .* tCycle;
end
SpectralSeparation = sum(sumSpectralSeparation) / (info.subSwath.numOfBursts *info.subSwath.samplesPerBurst);
end
function [ dopplerRate] = computeDopplerRate( info,burstIndex )
RangeDependentDopplerRate=computeRangeDependentDopplerRate(info.subSwath,burstIndex);
azTime = (info.sensor_para.first_line_time + info.sensor_para.last_line_time)/2.0;
velocity=getPositionVelocity(azTime,info);
v = getVelocity(velocity); 
steeringRate = info.subSwath.azimuthSteeringRate * (pi / 180.0);
krot = 2*v*steeringRate/info.sensor_para.Wavelength;
dopplerRate= RangeDependentDopplerRate * krot...
                            ./ (RangeDependentDopplerRate - krot);
end
function RangeDependentDopplerRate = computeRangeDependentDopplerRate(subSwath,burstIndex)
azFmRateList = subSwath.azFmRateList;
x = 1:subSwath.samplesPerBurst;
slrt = getSlantRangeTime(x,subSwath)*2;
dt = slrt - azFmRateList(burstIndex).t0;
RangeDependentDopplerRate =...
   azFmRateList(burstIndex).c0 + azFmRateList(burstIndex).c1.*dt + azFmRateList(burstIndex).c2.*dt.^2;
end
function slrt=getSlantRangeTime(x,subSwath)
 slrt=subSwath.slrTimeToFirstPixel +(x-1) * subSwath.rangePixelSpacing / 299792458;
end
function [Velocity]=getVelocity(POS)
Velocity=sqrt(POS.x*POS.x+POS.y*POS.y+POS.z*POS.z);
end
function velocity=getPositionVelocity(time,info)
nv = 8;
orbitStateVectors=info.orbit_para.orbit_points;
[m,~]=size(orbitStateVectors);
dt = (orbitStateVectors(m,1)-orbitStateVectors(1,1)) / (m-1);
            if (m <= nv) 
                i0 = 0;
                iN = m - 1;
            else 
                i0 = max(floor(((time - orbitStateVectors(1,1) )/ dt)) - nv/2 + 1, 0);
                iN = min(i0 + nv - 1, m - 1);
                if iN < m - 1
                i0 = i0;
                else
                i0=iN - nv + 1;
                end
            end
                position.x = 0;
                position.y = 0;
                position.z = 0;
                velocity.x = 0;
                velocity.y = 0;
                velocity.z = 0;
            for i = i0:iN
                orbI.time_mjd= orbitStateVectors(i+1,1);
                orbI.x_pos = orbitStateVectors(i+1,2);
                orbI.y_pos = orbitStateVectors(i+1,3);
                orbI.z_pos = orbitStateVectors(i+1,4);
                orbI.x_vel = orbitStateVectors(i+1,5);
                orbI.y_vel = orbitStateVectors(i+1,6);
                orbI.z_vel = orbitStateVectors(i+1,7);
                weight = 1;
                for j = i0:iN
                    if (j ~= i) 
                        time2 = orbitStateVectors(j+1,1);
                        weight =weight * (time - time2) / (orbI.time_mjd - time2);
                    end
                end
                    position.x = position.x+ weight * orbI.x_pos;
                    position.y = position.y+ weight * orbI.y_pos;
                    position.z = position.z +weight * orbI.z_pos;
                    velocity.x = velocity.x + weight * orbI.x_vel;
                    velocity.y = velocity.y +weight * orbI.y_vel;
                    velocity.z = velocity.z + weight * orbI.z_vel;
            end
end
function  [order] = sequential_netwrok(stacksize,N)
%N: degree
%stack size
m=1;
for i=1:stacksize-1
    for j=i+1:i+N
        if j>stacksize
            continue;
        end
   order(m,1)=i;order(m,2)=j; 
   m=m+1;
    end
end
end
function ifg = earlymultilooking(bm,bs,winsize)


if nargin < 3
    winsize = [7 7];
end

ifg = bm.*conj(bs);

h = ones(winsize(1),winsize(2));
ifg = filter2(h,ifg);
ifg = exp(1j*angle(ifg));
 
%     ifg1 = earlymultilooking(b1m,b1s);
% ifg2 = earlymultilooking(b2m,b2s);
% ifg = ifg1.*conj(ifg2);
% assuming that coh estimated by 2K sample with [15 15 ] winsize, then
% ifg = ifg.*coh.^2;
% eta = sum(ifg(a:b,a:b))/sum(coh(a:b,a:b).^2)
end
function [misrad,maxcoh] = periodogram(diffph) 
tempph = single(-1.5:0.02:1.5);
coh = zeros(length(tempph),1,'single');
for i = 1:length(tempph)
coh(i,1) = abs(real(mean(exp(1j*(diffph-tempph(i))))));
end
ind = find(coh == max(coh));
maxcoh = max(coh);
misrad = tempph(ind(1));
end