Add support of L4/5 fft filter

[alex-s-gardner]

L4 & 5 imagery is plagued by banding issues that are especially noticeable over low contrast snow and ice surfaces. This causes the feature tracking to lock onto banding discontinuities instead of surface features. To reduce the influence of banding we need to implement an fft filter that deconstructs the image, removes power in the across-track image direction then reforms the image.

Here is a Matlab that removes the banding:

function Ix = arImgFFT(validDomain, Ix)

  % remove across-track stripes in image

  % find the 4 corners of the image

  [m, n] = size(validDomain);
  s = regionprops(validDomain,'centroid');

  % sometimes there is more than one region (isloated
  % valid pixels) so pick largest region
  if length(s) > 1
      a = regionprops(validDomain,'area');
      idx = max([a.Area]) == [a.Area];
      s = s(idx);
      clear a idx
  end
  foo = false([m,n]);
  foo(round(s.Centroid(2)),round(s.Centroid(1))) = true;
  foo = bwdist(foo);
  foo(~validDomain) = 0;

  mC = round(m/2);
  nC = round(n/2);

  % bottom left
  [blR, blC] = find((max(max(foo(1:mC,1:nC))) == foo(1:mC,1:nC)),1, 'first');

  % bottom right
  [brR, brC] = find((max(max(foo(1:mC,(nC+1):end))) == foo(1:mC,(nC+1):end)),1, 'first');
  brC = brC+nC;

  % top left
  [tlR, tlC] = find((max(max(foo((mC+1):end,1:nC))) == foo((mC+1):end,1:nC)),1, 'first');
  tlR = tlR + mC;

  % top right
  [trR, trC] = find((max(max(foo((mC+1):end,(nC+1):end))) == foo((mC+1):end,(nC+1):end)),1, 'first');
  trR = trR + mC;
  trC = trC + nC;

  % figure
  % imagesc(validDomain)
  % hold on
  % plot(blC,blR, 'r*');
  % plot(brC,brR, 'g*');
  % plot(tlC, tlR, 'b*');
  % plot(trC,trR,  'k*');
  % plot(s.Centroid(1),s.Centroid(2),  'c*');

  % deteremine distance between corner points
  corners = [blC, blR; brC, brR; tlC, tlR; trC, trR];
  sep = pdist(corners,'euclidean');

  if any(sep < mC)
      % corners are likely in middle... try diferent
      % approach

      % bottom left
      [blR, blC] = find((max(max(foo(1:(round(mC/2)),1:end))) == foo(1:(round(mC/2)),1:end)), 1, 'first');

      % bottom right
      [brR, brC] = find((max(max(foo(1:end,round(nC/2*3):end))) == foo(1:end,round(nC/2*3):end)), 1, 'first');
      brC = brC+round(nC/2*3)-1;

      % top left
      [tlR, tlC] = find((max(max(foo(round(mC/2*3):end,1:end))) == foo(round(mC/2*3):end,1:end)),1, 'first');
      tlR = tlR + round(mC/2*3) -1;

      % top right
      [trR, trC] = find((max(max(foo(1:end,1:round(nC/2)))) == foo(1:end,1:round(nC/2))), 1, 'first');

      corners = [blC, blR; brC, brR; tlC, tlR; trC, trR];
      sep = pdist(corners,'euclidean');
      if any(sep < mC)
          % at least 2 points are too close to each other
          error('two or more recovered image corner locations are too close to eachother\n add to cleanLandsatDataDir corruptImages list then run cleanLandsatDataDir:\n %s'...
              , fooFileName)
      end

  end

  clear foo

  slope1 = atan((brR - blR) / (brC - blC));
  slope2 = atan((trR - tlR) / (trC - tlC));

  slope3 = atan((tlR - blR) / (blC - tlC));
  slope4 = atan((brR - trR) / (trC - brC));

  fooA = false(size(validDomain));
  center = round(size(fooA)/2);

  fooA(floor(center(1)) - 70 : floor(center(1)) + 70, :) = 1; % create mask - horizontal line
  fooA(:,(floor(center(2)) - 100) : (floor(center(2)) + 100)) = 0; % with zeros at the center

  % slope will be reduced is image is truncated at corners
  % so choose max slope
  A = imrotate(fooA,rad2deg(max([slope1,slope2])), 'crop');
  B = imrotate(fooA,rad2deg(max([slope3,slope4])), 'crop');

  % shift mask to center
  shiftCtr = round([s.Centroid(2), s.Centroid(1)]) - center ;
  tform = affine2d([1 0 0; 0 1 0; shiftCtr(2) shiftCtr(1) 1]);
  [A,~] = imwarp(A,tform);
  [B,~] = imwarp(B,tform);

  %Perform 2D FFTs
  for k = 1:length(size(Ix,3))

      fftIm = Ix(:,:,k);
      fftIm(fftIm>3) = 3;
      fftIm(fftIm<-3) = -3;
      fftIm(isnan(fftIm)) = 0;
      fftIm = fftshift(fft2(double(fftIm)));
      P = abs(fftIm);
      mP =  mean(P(:));
      stdP = std(P(:));
      P = (P-mP) > 10*stdP;
      sA = sum(P(A));
      sB = sum(P(B));

      % imagesc(P+A)

      if (sA/sB >= 2 || sB/sA >= 2) && (sA > 500 || sB > 500)
          if sA > sB
              mask = A;
          elseif sB > sA
              mask = B;
          end

          %figure; imshow(real(ifft2(ifftshift(fftIm .* (1-(mask ))))));
          foo1 = isnan(Ix(:,:,k));
          foo2 = real(ifft2(ifftshift(fftIm .* (1-(mask)))));
          foo2(foo1) = nan;

          %figure; imshow(foo2);caxis([-1.5 1.5])

          Ix(:,:,k) = foo2;

          clear fftIm mask foo1 foo2
      end
  end

[jhkennedy]

Notes:

• We'll be processing Band 2 typically for L4,5
• This filter would be applied after the Wallis filter is applied

[jhkennedy]

Done in #66