nasa-jpl / autoRIFT

A Python module of a fast and intelligent algorithm for finding the pixel displacement between two images
Apache License 2.0
212 stars 52 forks source link

Need to add additional filter to support Landsat 7 SLC-off imagery #52

Closed alex-s-gardner closed 2 years ago

alex-s-gardner commented 2 years ago

In the initial preprocessing of the imagery (e.g. application of the Wallis filter) we need to add an option for L7 filtering which requires filling the missing data stripes with random noise that as the same mean and std as the filtered valid L7 imagery. Below is the prototype code written in Matlab:

% add random noise to SLC_off imagery
% Landsat 7 data aquired AFTER 2003, 5, 31

%% ----------------- READ IN IMAGERY ---------------------
% check if image is ETM+ SLC-off
I = imread('data/LE07_L1TP_009011_20130809_20200907_02_T1_B8.tif');
I = single(I);

WallisFilterWidth = 5; % filter width <- already implemented in autoRIFT
StandardDeviationCutoff = 0.25; <- already implemented in autoRIFT

%% ----------------- START OF NEW FILTER ---------------------
% make random number generator unique
ImNoData = 0;
rng(round(sum(pwd)+now*100),'twister');

% buffer values lines by sqrt(((filtWidth-1)/2)^2) + 0.01 for
% rounding
buff = sqrt(2*((WallisFilterWidth-1)/2)^2) + 0.01;

validDomain = I  ~= ImNoData; % validDomain is the image frame
ImNoData = single(ImNoData);

% find edges of image, this makes missing scan lines valid and will
% later be filled with random white noise
missingData = bwdist(validDomain,'euclidean')<30; % missingData is missing data within image frame
missingData = missingData & ~validDomain;
missingData = bwdist(missingData,'euclidean')<= buff;

% trying to frame out the image
validDomain = validDomain | missingData;

LowStd = stdfilt(I, ones(WallisFilterWidth))<StandardDeviationCutoff;
LowStd = bwdist(LowStd,'euclidean')<= buff;
missingData(:,:,1) = (missingData(:,:,1) | LowStd) & validDomain;

%% ----------------- This is already coded in autoRIFT ----------------- 
% apply high pass filter
colFiltParallel = true;
fun = @(block_struct) wallis_filt(block_struct.data, WallisFilterWidth);
I = blockproc(I,[10000 10000], fun, ...
    'BorderSize', [WallisFilterWidth, WallisFilterWidth],'PadPartialBlocks', ...
    false, 'TrimBorder', true, 'UseParallel', colFiltParallel);

% if scan line corrector off image
% scale to fit range of int8

validData = validDomain & ~missingData;
I(~validData) = 0;

% populate random values with same std and mean
fillData = validDomain & missingData;
R = randn([sum(fillData(:)),1]);
I(fillData) = R;