r_interp_opc.m 10.3 KB
function r_interp_opc(expNumber,pretreat,zone,color,iFirst,nSkip,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)
% 
% Syntax:  r_interp_opc(expNumber,zone,color,iFirst,nSkip,iEnd)
% 
% Inputs:
%    expNumber - The date number (yyyymmdd). Points to a directory
%    pretreat  - 0: distorted images, 1: undistorted images, 2: stabilized,
%    3: undist + stable
%    zone      - The part of the scene to analyse. See the code for details
%    color     - 0: grayscale, 1: RGB image,
%    iFirst    - The first image to be treated
%    nSkip     - The number of images skipped
%    iEnd      - The last image to be treated
%                (if iEnd is set to zero, the loop is done from iFirst to
%                 nFiles)
% 
% Example:
%    r_interp_opc(20140203,1,3,0,1,1,0)
% 
% Other m-files required: g_rect toolbox
%   (http://demeter.uqar.ca/g_rect)
% Subfunctions: none
% MAT-files required: the rectification parameter file
% 
% See also: G_RECT, R_UNDISTORT, R_ROTATE

% Author: Dany Dumont
% UQAR/ISMER
% email: dany_dumont@uqar.ca
% Website: http://www.ismer.ca/dumont-dany
% February 2014
% ______________________________________________________________________

estr = num2str(expNumber);
zstr = num2str(zone);

rootdir = '/home/nicopa01/data/opc/images';
expdir  = [rootdir,'/',estr];
if pretreat == 0;
    origdir = [expdir,'/original'];
    intdir  = [expdir,'/interp'];
    calfile = [rootdir,'/calib/opc_rect_' estr '.mat']; 
elseif pretreat == 1;
    origdir = [expdir,'/original_undist'];
    intdir  = [expdir,'/interp_undist'];
    calfile = [rootdir,'/calib/opc_rect_' estr '_undist.mat']; 
elseif pretreat == 2;
    origdir = [expdir,'/original_stable'];
    intdir  = [expdir,'/interp_stable'];
    calfile = [rootdir,'/calib/opc_rect_' estr '_stable.mat']; 
elseif pretreat == 3;
    origdir = [expdir,'/undist_stable'];
    intdir  = [expdir,'/interp_undist_stable'];
    calfile = [rootdir,'/calib/opc_rect_' estr '_undist_stable.mat']; 
end
zonedir = [intdir,'/zone',zstr];
tifdir  = [zonedir,'/tif'];  

%% Load the calibration file
load(calfile)

%% Creation of the grid


if zone == 1            % Zone 1 (à refaire)
    dx = 0.000009;      % degree_lon ~  7 m
    dy = 0.000009;      % degree_lat ~ 11 m
    lon1 = -68.824478; lon2 = -68.825343;
    lat1 = 48.383910; lat2 = 48.398269;
elseif zone == 2        % Zone 2 (à refaire)
    dx = 0.00002;       % degree_lon ~ 1.4 m
    dy = 0.00002;       % degree_lat ~ 2.2 m
    lon1 = -68.8650; lon2 = -68.8250;
    lat1 =  48.3340; lat2 =  48.3660;
elseif zone == 3        % Zone 3 (à refaire)
    dx = 0.000002;      % degree_lon ~ 1.4 m
    dy = 0.000002;      % degree_lat ~ 2.2 m
    lon1 = -68.8400; lon2 = -68.8310;
    lat1 =  48.3385; lat2 =  48.3430;
elseif zone == 4        % Zone 4 (entrée de la Baie du Ha Ha!)
    dx = 0.000006;      % degree_lon ~ 1.4 m
    dy = 0.000006;      % degree_lat ~ 2.2 m
    lon1 = -68.8460; lon2 = -68.8310;
    lat1 =  48.3385; lat2 =  48.3510;
elseif zone == 5        % Zone 5 (à refaire)
    dx = 0.000006;      % degree_lon ~ 1.4 m
    dy = 0.000006;      % degree_lat ~ 2.2 m
    lon1 = -68.8460; lon2 = -68.8250;
    lat1 =  48.3385; lat2 =  48.3510;
elseif zone == 6        % Zone 6 (à refaire)
    dx = 0.000003;      % degree_lon ~  m
    dy = 0.000004;      % degree_lat ~  m
    lon1 = -68.8400; lon2 = -68.8320;
    lat1 =  48.3385; lat2 =  48.3490;
elseif zone == 7        % Zone 7 (entrée baie + au large)
    dx = 0.000009;      % degree_lon ~  m
    dy = 0.000012;      % degree_lat ~  m
    lon1 = -68.862252; lon2 = -68.825343;
    lat1 = 48.338220; lat2 = 48.377713 ;
elseif zone == 8        % Zone 8 (entrée baie du Ha Ha!)
    dx = 0.000009;      % degree_lon ~  m
    dy = 0.000009;      % degree_lat ~  m 0.000009
    lon1 = -68.862252; lon2 = -68.825343;
    lat1 = 48.338220; lat2 = 48.356827;
elseif zone == 9        % Zone 9 (Chenal du Bic)
    dx = 0.00002;       % degree_lon ~  m
    dy = 0.00002;       % degree_lat ~  m 
    lon1 = -68.89; lon2 = -68.81;
    lat1 = 48.336168; lat2 = 48.40;
elseif zone == 10        % Zone 10 (autour de l'awac 2015)
    dx = 0.00004;       % degree_lon ~  m
    dy = 0.00004;       % degree_lat ~  m 
    lon1 = -68.8850; lon2 = -68.8700;
    lat1 = 48.34000; lat2 = 48.3550;
end

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);

cd(origdir)

imfiles = dir('*.jpg');
nFiles  = size(imfiles,1);

%% Masking
mask = ones(5472,3648);

if zone == 1
    % Zone 1 (à refaire)
    mask(:,   1:1550) = NaN;         % above the horizon
    mask(:,2500:end ) = NaN;         % below the tree line
elseif zone == 2
    % Zone 2 (à refaire)
    mask(:,   1:1550) = NaN;         % above the center of the channel
    mask(:,2500:end ) = NaN;         % below the tree line
    mask(1:1500,:)    = NaN;
elseif zone == 3
    % Zone 3 (à refaire)
    mask(:,   1:1750) = NaN;         % above the center of the channel
    mask(:,2500:end ) = NaN;         % below the tree line
    mask(1:2500,:)    = NaN;
elseif zone == 4 || zone == 5 || zone == 6
    % Zone 4 to 6 (entrée de la Baie du Ha Ha!)
    mask(:,   1:1550) = NaN;         % above the center of the channel
    mask(:,2500:end ) = NaN;         % below the tree line
    mask(1:2300,:)    = NaN;
elseif zone == 7 || zone == 8
    % Zone 7 & 8 (Baie du Ha Ha! et large)
    mask(:,   1:1660) = NaN;         % above the center of the channel
    mask(:,2500:end) = NaN;         % below the tree line
    mask([1:1000 5000:end],:)    = NaN;
elseif zone == 9
    % Zone 9 (chenal du Bic)
    mask(:,   1:1660) = NaN;         % above the center of the channel
    mask(:,2500:end) = NaN;         % below the tree line
    mask([1:1000 5000:end],:)    = NaN;
elseif zone == 10
    % Zone 10 (autour de l'AWAC du Bic 2015)
    mask(:,   1:1660) = NaN;         % above the center of the channel
    mask(:,2500:end) = NaN;         % below the tree line
    mask([1:1000 5000:end],:)     = NaN;
else
    
    disp( ' Error :  Choose a valid zone number');
    return
end

%%

LON = mask.*LON; %#ok<NODEF>
LAT = mask.*LAT; %#ok<NODEF>

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

if ~exist(intdir,'dir')
    mkdir(intdir)
end

if ~exist(zonedir,'dir')
    mkdir(zonedir)
end

if ~exist(tifdir,'dir')
    mkdir(tifdir)
end

save([zonedir,'/zone',zstr,'.mat'],'X1','Y1');

if iEnd == 0
    iEnd = nFiles;
end

for n = iFirst:nSkip:iEnd
    
    imfilename = imfiles(n).name;
    imNumber   = imfilename(1:13);

    I0 = imread(imfilename);
    I0 = double(I0);
    
     if color == 1
        
        R0 = mask.*I0(:,:,1)';
        G0 = mask.*I0(:,:,2)';
        B0 = mask.*I0(:,:,3)';
        
        r0 = R0(:);
        g0 = G0(:);
        b0 = B0(:);
        
        r0(isnan(r0)) = [];
        g0(isnan(g0)) = [];
        b0(isnan(b0)) = [];
        x1(isnan(x1)) = [];
        y1(isnan(y1)) = [];
        
        Fr = TriScatteredInterp(x1,y1,r0);
        Fg = TriScatteredInterp(x1,y1,g0);
        Fb = TriScatteredInterp(x1,y1,b0);
        
        R = flipud(Fr(X1,Y1));
        G = flipud(Fg(X1,Y1));
        B = flipud(Fb(X1,Y1));
        
        I(:,:,1) = R;
        I(:,:,2) = G;
        I(:,:,3) = B;
        
        I = uint8(I);

        if pretreat == 0;
            filename = [zonedir,'/',imNumber,'_',zstr,'_rbg.jpg'];
        elseif pretreat == 1;
            filename = [zonedir,'/',imNumber,'_',zstr,'_rbg_undist.jpg'];
        elseif pretreat == 2;
            filename = [zonedir,'/',imNumber,'_',zstr,'_rbg_stable.jpg'];
        elseif pretreat == 3;
            filename = [zonedir,'/',imNumber,'_',zstr,'_undist_stable.jpg'];
        end
        
        imwrite(I,filename,'jpg');
        
        % create geotiff file
        R = georasterref();
        R.RasterSize = [size(I)];
        R.Latlim = [lat1 lat2];
        R.Lonlim = [lon1 lon2];
        R.ColumnsStartFrom = 'north';
        
        if pretreat == 1;
            filename = [zonedir,'/tif/',imNumber,'_',zstr,'_rbg_undist.tif'];
        elseif pretreat == 0;
            filename = [zonedir,'/tif/',imNumber,'_',zstr,'_rbg.tif'];
        elseif pretreat == 2;
            filename = [zonedir,'/tif/',imNumber,'_',zstr,'_rbg_stable.tif'];
        elseif pretreat == 3;
            filename = [zonedir,'/tif/',imNumber,'_',zstr,'_undist_stable.tif'];
        end
        
        geotiffwrite(filename, I, R);
        
    elseif color == 0
        
        if pretreat == 2
            K0 = mask.*I0';
        else
            K0 = mask.*mean(I0,3)';
        end
        
        k0 = K0(:);
        
        k0(isnan(k0)) = [];
        x1(isnan(x1)) = [];
        y1(isnan(y1)) = [];
        
        Fk = TriScatteredInterp(x1,y1,k0);
        
        K = flipud(Fk(X1,Y1));
        
        I = uint8(K);
        
        if pretreat == 0;
            filename = [zonedir,'/',imNumber,'_',zstr,'_bw.jpg'];
        elseif pretreat == 1;
            filename = [zonedir,'/',imNumber,'_',zstr,'_bw_undist.jpg'];
        elseif pretreat == 2;
            filename = [zonedir,'/',imNumber,'_',zstr,'_bw_stable.jpg'];
        elseif pretreat == 3;
            filename = [zonedir,'/',imNumber,'_',zstr,'_undist_stable.jpg'];
        end
        
        
        imwrite(I,filename,'jpg');
        
        % create geotiff file
        R = georasterref();
        R.RasterSize = [size(I)];
        R.Latlim = [lat1 lat2];
        R.Lonlim = [lon1 lon2];
        R.ColumnsStartFrom = 'north';
        
        if pretreat == 0;
            filename = [zonedir,'/tif/',imNumber,'_',zstr,'_bw.tif'];
        elseif pretreat == 1;
            filename = [zonedir,'/tif/',imNumber,'_',zstr,'_bw_undist.tif'];
        elseif pretreat == 2;
            filename = [zonedir,'/tif/',imNumber,'_',zstr,'_bw_stable.tif'];
        elseif pretreat == 3;
            filename = [zonedir,'/tif/',imNumber,'_',zstr,'_undist_stable.tif']; 
        end
        
        geotiffwrite(filename, I, R);
    end
    
        
    disp([' ',imNumber]);
    
end