%
% Purpose: Fourier-based image destriping
% Written by: H.E.
% Date: 2/2009
%

close all;
clear;
EV_250m = double(hdfread('MOD01.A2009035.1115.hdf', '/EV_250m', 'Index', {[1  1],[1  1],[1716  1248]}));
figure; imagesc(EV_250m); colormap('gray');
title('Terra MODIS L1A (from day Feb. 4, 2009) 645 \mum or 856 \mum');

im_fft=fftshift(fft2(EV_250m));
figure; imagesc(log(abs(im_fft))); colormap(gray);
[x,y]=size(EV_250m); xx=1:x;
ln=log(im_fft(:,y/2+1)); p=polyfit(xx,ln',2); ln_e=1.75+polyval(p,xx);
idx=find(ln'>ln_e); idx=[idx(1:31) idx(63:length(idx))];

figure; plot(xx,ln); grid; hold; plot(xx,ln_e); hold;
figure; ln(idx)=0.0; plot(xx,ln);

for i=2:length(idx)-1
    im_fft(idx(i)-1:idx(i)+1,y/2:y/2+2)=0.0;
end
%figure; imagesc(log(abs(im_fft))); colormap(gray);

im_f=abs(ifft2(fftshift(im_fft)));
imagesc(im_f); colormap(gray);
title('Terra MODIS L1A Destriped 645 \mum or 856 \mum');