findMatdetecFootprint.asv

function Matdetec2 = findMatdetecFootprint(DetectFootPrint,XPsizeOut,YPsizeOut)
% Fixed the calculation of footprint of Sentinel-2 detector based on old
% version of Sentinel-2 data (17/12/2021 Shi)

clear Matdetec k
% sorted foot print from 1 to 4
[DetectFootPrint.Ndect, idnews] = sort(DetectFootPrint.Ndect);
DetectFootPrint.Nband = DetectFootPrint.Nband(idnews);
DetectFootPrint.Data = DetectFootPrint.Data(idnews);

for i = 1:length(DetectFootPrint.Ndect)
  
  IDX = knnsearch(XPsizeOut',DetectFootPrint.Data{i,1}(:,1));
  IDY = knnsearch(YPsizeOut,DetectFootPrint.Data{i,1}(:,2));
  
  dum2 = single(poly2mask(double(IDX), double(IDY),length(XPsizeOut),length(XPsizeOut))) ;
  clear IDX IDY;
  dum2 = conv2(dum2,ones(3),'same')>0; % to fill boundary
  Matdetec(:,:,i) = dum2;
  clear dum*
  
  %   find orientation of detect + slope for computing perpendicular kernel
  I=nan(size(Matdetec,2),1);
  for ii=1:size(Matdetec,1)
    dum=find(Matdetec(ii,:,i)==1,1,'first');
    if ~isempty(dum)
      I(ii,1)=dum;
    end
  end
  clear dum;
  J = [1:size(Matdetec,1)]' ;
  test = ~isnan(I) & I > 1 & I < size(Matdetec,2);
  warning off all % if warning => not enough point => slope=0 => good because tile boundary
  k{i,1} = polyfit(J(test),I(test),1);
  clear test;
    
  I=nan(size(Matdetec,2),1);
  for ii=1:size(Matdetec,1)
    dum=find(Matdetec(ii,:,i)==1,1,'last');
    if ~isempty(dum)
      I(ii,1)=dum;
    end
  end
  J = [1:size(Matdetec,1)]' ;
  test = ~isnan(I) & I > 1 & I < size(Matdetec,2);
  k{i,2} = polyfit(J(test),I(test),1);
  clear test;
  warning on all

end

% mediane
for i = 1:length(DetectFootPrint.Ndect)-1
  mediane = mean( [k{i,2} ; k{i+1,1}] ) ;
  
  k{i,2} = mediane ;
  k{i+1,1} = mediane ;
  clear mediane;
end
J = [1:size(Matdetec,1)]' ;
I = [1:size(Matdetec,2)] ;

[Jmat Imat] = meshgrid(I,J);
clear I J;

Matdetec2 = nan(size(Matdetec,1),size(Matdetec,2));
clear Matdetec;
for i = 1:length(DetectFootPrint.Ndect)
  
  liminf = polyval(k{i,1},Jmat);
  limsup = polyval(k{i,2},Jmat);
  
  if sum(k{i,2}) == 0 % footprint at low-right corner
      Matdetec2(Imat>=liminf) = DetectFootPrint.Ndect(i) ;
  elseif sum(k{i,1}) == 0 % footprint at up-left corner
     Matdetec2(Imat<=limsup) = DetectFootPrint.Ndect(i) ;
  else
     Matdetec2(Imat>=liminf & Imat<=limsup) = DetectFootPrint.Ndect(i) ;
  end
  clear liminf limsup;
end
clear Imat ImatJ k;

Matdetec2 = Matdetec2';