Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@jasonleepolyu
jasonleepolyu committed Jun 5, 2017
1 parent 61805ae commit 718100d
Show file tree
Hide file tree
Showing 2 changed files with 225 additions and 0 deletions.
52 changes: 52 additions & 0 deletions ProPADS15to20.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,52 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% bilinear and cubic convolution resampling %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all; close all; clc

%% Input original image
%%%% The image 'PXS_15m_B1_B7_Brovey_odd.tif' was generated by implementing 'ProPADS30to15.m'

[I0, R] = geotiffread('PXS_15m_B1_B7_Brovey_odd.tif');

info = geotiffinfo('20m_B4.tif'); % Sentinel-2 20 m reference band, just its 'info' and 'R' were used here.
[temp, R] = geotiffread('20m_B4.tif');

b1 = I0(:,:,1);
b2 = I0(:,:,2);
b3 = I0(:,:,3);
b4 = I0(:,:,4);
b5 = I0(:,:,5);
b6 = I0(:,:,6);
b7 = I0(:,:,7);
%% Input scale factor
K = 15/20; % projection of Landsat 15 m data to Sentinel-2 20 m resolution
% K = 30/20; % projection of Landsat 30 m data to Sentinel-2 20 m resolution
%% cubic convolution (cc) or bilinear (bl) resampling
resampler = 'cc'; % or resampling_type = 'bl';

b1 = bl_cc_resampling_at(b1,K,resampler);
b2 = bl_cc_resampling_at(b2,K,resampler);
b3 = bl_cc_resampling_at(b3,K,resampler);
b4 = bl_cc_resampling_at(b4,K,resampler);
b5 = bl_cc_resampling_at(b5,K,resampler);
b6 = bl_cc_resampling_at(b6,K,resampler);
b7 = bl_cc_resampling_at(b7,K,resampler);

%% Output
Output(:,:,1) = b1;
Output(:,:,2) = b2;
Output(:,:,3) = b3;
Output(:,:,4) = b4;
Output(:,:,5) = b5;
Output(:,:,6) = b6;
Output(:,:,7) = b7;
[row, col, band] = size(b5)
R.RasterSize = [row col];
% geotiffwrite('registered_20m.tif',int16(Output),R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite('registered_20m.tif',int16(Output),R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

% Output = uint8(Output(:,:,[7,5,4]));
% figure, imshow(Output(:,:,[7,5,4]));



173 changes: 173 additions & 0 deletions bl_cc_resampling_at.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,173 @@
function Output = bl_cc_resampling_at(I, scale_factor,resampler)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% I: the original image or bands
%%% scale_factor
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% pad image by duplicate the last row and column
[nrow, ncol] = size(I);
I = [I; I(nrow, :); I(nrow, :); I(nrow, :);];
I = [I I(:, ncol) I(:, ncol) I(:, ncol)];

%% Initialize the output image
nrow_new = nrow*scale_factor;
ncol_new = ncol*scale_factor;
Output = zeros(nrow_new, ncol_new);

%% Ratio of output image to the original image
% y_ratio = (nrow)/(nrow*scale_factor);
% x_ratio = (ncol)/(ncol*scale_factor);

%% BL and CC Affine transformation coefficients
%% landslide
% a0 = -1.346995354 + 12.5/20; % +12.5/20 col shift
% a1 = 1.000008878;
% a2 = 0.000033128;
% b0 = -1.218221657 + 12.5/20; % +12.5/20 row shift
% b1 = -0.000042180;
% b2 = 1.000012239;
%
% col_shift = 12.5; % geolocation shift (upper left) between Landsat 15 m and Sentinel-2 20 m data
% row_shift = 12.5;

%% New_york

% a0 = -0.771690110 + 12.5/20; % +12.5/20 col shift
% a1 = 0.999915075;
% a2 = -0.000189602;
% b0 = 1.078519619 + 2.5/20; % +2.5/20 row shift
% b1 = -0.000064901;
% b2 = 0.999881731;

% col_shift = 12.5; % geolocation shift (upper left) between Landsat 15 m and Sentinel-2 20 m data
% row_shift = 2.5;

%% crop field

% a0 = -1.007142445 + 12.5/20; % +12.5/20 col shift
% a1 = 1.000007977;
% a2 = 0.000012005;
% b0 = 0.225335187 + 2.5/20; % +2.5/20 row shift
% b1 = 0.000015280;
% b2 = 0.999996748;

% col_shift = 12.5; % geolocation shift (upper left) between Landsat 15 m and Sentinel-2 20 m data
% row_shift = 2.5;

%% burned_area

% a0 = -0.772688990 + 2.5/20; % +2.5/20 col shift
% a1 = 1.000018270;
% a2 = 0.000036507;
% b0 = 0.489724920 + 2.5/20; % +2.5/20 row shift
% b1 = -0.000030781;
% b2 = 1.000013579;

% col_shift = 2.5; % geolocation shift (upper left) between Landsat 15 m and Sentinel-2 20 m data
% row_shift = 2.5;

%% forest

a0 = 0.184715799 + 2.5/20; % col shift
a1 = 0.999815008;
a2 = -0.000100652;
b0 = 0.221510137 + 2.5/20; % row shift
b1 = -0.000149049;
b2 = 0.999983763;

col_shift = 2.5; % geolocation shift (upper left) between Landsat 15 m and Sentinel-2 20 m data
row_shift = 2.5;

for j = 1:nrow_new

for i = 1:ncol_new

m = a0 + a1 * i + a2 * j; % col affine transformation
n = b0 + b1 * i + b2 * j; % row affine transformation

if (m >=1 && m < ncol_new && n >=1 && n < nrow_new )
%% without considering the geolocation shift between Landsat and Sentinel-2 data
% YY = floor((n-1)*y_ratio+1); % row
% y = (n-1)*y_ratio+1 - YY;
%
% XX = floor((m-1)*x_ratio+1); % col
% x = (m-1)*x_ratio+1 - XX;
%% considering the geolocation shift between Landsat and Sentinel-2 data
floaty = ((n-0.5)*20 - row_shift)/15+1; % 15 denotes 15 m to 20 m; 30 can be used for 30 m to 20 m;
if (floaty - floor(floaty))>=0.5
YY = floor(floaty); % row
y = floaty - YY - 0.5;
else
YY = floor(floaty)-1; % row
y = floaty - YY - 0.5 ;
end

floatx = ((m-0.5)*20 - col_shift)/15+1; % 15 denotes 15 m to 20 m; 30 can be used for 30 m to 20 m;
if (floatx - floor(floatx))>=0.5
XX = floor(floatx); % col
x = floatx - XX - 0.5;
else
XX = floor(floatx)-1; % col
x = floatx - XX - 0.5;
end
if y<0 || x<0
fprintf('shit this is wrong ');
end


if (YY >1 && YY < nrow && XX >1 && XX < ncol )
%% CC and BL resampling
switch resampler
case 'cc'
% cubic convolution

A = [cc_kernel(y + 1) cc_kernel(y + 0) cc_kernel(y - 1) cc_kernel(y - 2)];

B = [I(YY-1, XX-1) I(YY-1, XX+0) I(YY-1, XX+1) I(YY-1, XX+2);
I(YY+0, XX-1) I(YY+0, XX+0) I(YY+0, XX+1) I(YY+0, XX+2);
I(YY+1, XX-1) I(YY+1, XX+0) I(YY+1, XX+1) I(YY+1, XX+2);
I(YY+2, XX-1) I(YY+2, XX+0) I(YY+2, XX+1) I(YY+2, XX+2);];

C = [cc_kernel(x + 1) cc_kernel(x + 0) cc_kernel(x - 1) cc_kernel(x - 2)]';

Output(j,i) = double(A)*double(B)*double(C);

case 'bl'
% bilinear

Output(j,i) = (I(YY+0, XX+0)*(1-x)+I(YY+0, XX+1)*x)*(1-y)+ (I(YY+1, XX+0)*(1-x)+I(YY+1, XX+1)*x)*y;
end
end
end
end
end

end

function cc = cc_kernel(x)

% Equation (15) in Keys's paper "Cubic Convolution Interpolation for Digital Image Processing, 1981"

if abs(x)>=0 && abs(x)< 1
cc = 1.5 * abs(x)^3 - 2.5 * abs(x)^2 + 1;
end

if abs(x)>=1 && abs(x)< 2
cc = -0.5 * abs(x)^3 + 2.5 * abs(x)^2 - 4 * abs(x) + 2;
end

if abs(x)>=2
cc = 0;
end


end


% Output(j,i) = I(YY+0, XX+0);
% YY = floor(floaty); % row
% y = floaty - YY - 0.5;
% YY = floor(((n)*20-12.5)/15); % row
% y = ((n)*20-12.5)/15 - YY;
%
% XX = floor(((m)*20-12.5)/15); % col
% x = ((m)*20-12.5)/15 - XX;

0 comments on commit 718100d

Please sign in to comment.