r_interp_pitle_2016.m 2.74 KB
function r_interp_pitle_2016(cam,iBegin,iSkip,iEnd)
% R_INTERP_OPC - Uses the rectification data to interpolate RGB values on a
% regular grid and save it to an image. It is specifically adapted to
% images taken from the Observatoire du Pic Champlain (OPC)

rootdir = '/share/archives/dumoda01/data/photogrammetry/pitle/2016/';

%% image repertory

imfiles = dir([rootdir,cam,'/original/2016*']);

%% area masking
if cam == 'heroE'
	lon1 = -77.4398; lon2 = -76.2889;
	lat1 = 72.7179; lat2 = 72.9484;
	mask = ones(2976,1984);
	mask(:,1:800) = NaN;
	mask(:,1550:end) = NaN;

elseif cam ==  'heroW'
	lon1 = -79.3746; lon2 = -77.5490;
	lat1 = 72.5474; lat2 = 72.8927;
	mask = ones(6000,4000);
	mask(:,1:2000) = NaN;

elseif cam == 'buttP'
	lon1 = -76.4498; lon2 = -75.5091;
	lat1 = 72.5774; lat2 = 72.85739;
	mask = ones(5184,3456);
	mask(:,1:1500) = NaN;
	mask(:,2300:end) = NaN;

elseif cam == 'erikH'
   	lon1 = -76.8429; lon2 = -75.7753;
	lat1 = 72.5978; lat2 = 72.8511;
	mask = ones(5184,3456);
	mask(:,1:1500) = NaN;
	mask(:,2300:end) = NaN;

elseif cam == 'guysB'
	lon1 = -76.8429; lon2 = -75.7753;
	lat1 = 72.5978; lat2 = 72.8511;
	mask = ones(5517,3979);
	mask(:,1:1500) = NaN;
	mask(:,2300:end) = NaN;
    imfiles = dir([rootdir,cam,'/original/2016*_rotated.jpg']);
end

%% x/y resolution
dx = 0.0008; dy = 0.0002;
x1vec =  lon1:dx:lon2;
y1vec =  lat1:dy:lat2;
y1vec = y1vec';
xpolyvec = [lon1 lon1 lon2 lon2];
ypolyvec = [lat1 lat2 lat2 lat1];

[X1,Y1] = meshgrid(x1vec,y1vec);


%% data loading
load([ rootdir cam '/work/g_rect_' cam '_20160630.mat']);

LON = mask.*LON;
LAT = mask.*LAT;

x1 = LON(:); x2 = LON(:);
y1 = LAT(:); 

%% time retreaving
for i = 1:length(imfiles)
    imfiles(i).N = datenum(imfiles(i).name(1:15),'yyyymmdd_HHMMSS');
end

nBegin = datenum(num2str(iBegin),'yyyymmddHHMM');
nEnd   = datenum(num2str(iEnd),'yyyymmddHHMM');

N = [imfiles.N];
N = find(N>nBegin & N<nEnd);


%% image processing loop
for n = N(1):iSkip:N(end)
    imfilename = imfiles(n).name;
    imNumber   = imfilename(1:15);

    I0 = imread([ rootdir cam '/original/' imfilename ]);
    I0 = double(I0);
    
    i0 = mean(I0,3);
    K0 = mask.*i0'; 
    
    k0 = K0(:); 
    
    k0(isnan(x2)) = [];
    x1(isnan(x1)) = [];
    y1(isnan(y1)) = [];
   
    Fk = TriScatteredInterp(x1,y1,k0);
    
    K = flipud(Fk(X1,Y1));
    
    I = uint8(K);
    
    jpgname = [ rootdir cam '/work/jpg/' cam '_' imNumber '_interp.jpg'];
    tifname = [ rootdir cam '/work/tif/' cam '_' imNumber '.tif'];
        
	imwrite(I,jpgname,'jpg');
        
    % create geotiff file
    R = georasterref();
    R.RasterSize = [size(I)];
    R.Latlim = [lat1 lat2];
    R.Lonlim = [lon1 lon2];
    R.ColumnsStartFrom = 'north';
    
    geotiffwrite(tifname, I, R);
%    end
    
        
    disp([' ',imNumber]);
end