Blame view

r_interp_pitle_2016.m 2.74 KB
79657b4e   Paul Nicot   prise en compte d...
1
function r_interp_pitle_2016(cam,iBegin,iSkip,iEnd)
a48ac2ca   Paul Nicot   mfile pour le tra...
2
3
4
5
% 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)

5872298f   Paul Nicot   maj pour fonction...
6
rootdir = '/share/archives/dumoda01/data/photogrammetry/pitle/2016/';
a48ac2ca   Paul Nicot   mfile pour le tra...
7

79657b4e   Paul Nicot   prise en compte d...
8
9
10
11
12
%% image repertory

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

%% area masking
5872298f   Paul Nicot   maj pour fonction...
13
if cam == 'heroE'
a48ac2ca   Paul Nicot   mfile pour le tra...
14
15
16
17
18
19
	lon1 = -77.4398; lon2 = -76.2889;
	lat1 = 72.7179; lat2 = 72.9484;
	mask = ones(2976,1984);
	mask(:,1:800) = NaN;
	mask(:,1550:end) = NaN;

5872298f   Paul Nicot   maj pour fonction...
20
elseif cam ==  'heroW'
a48ac2ca   Paul Nicot   mfile pour le tra...
21
22
23
24
25
	lon1 = -79.3746; lon2 = -77.5490;
	lat1 = 72.5474; lat2 = 72.8927;
	mask = ones(6000,4000);
	mask(:,1:2000) = NaN;

5872298f   Paul Nicot   maj pour fonction...
26
elseif cam == 'buttP'
a48ac2ca   Paul Nicot   mfile pour le tra...
27
28
29
30
31
32
	lon1 = -76.4498; lon2 = -75.5091;
	lat1 = 72.5774; lat2 = 72.85739;
	mask = ones(5184,3456);
	mask(:,1:1500) = NaN;
	mask(:,2300:end) = NaN;

5872298f   Paul Nicot   maj pour fonction...
33
elseif cam == 'erikH'
a48ac2ca   Paul Nicot   mfile pour le tra...
34
35
36
37
38
39
   	lon1 = -76.8429; lon2 = -75.7753;
	lat1 = 72.5978; lat2 = 72.8511;
	mask = ones(5184,3456);
	mask(:,1:1500) = NaN;
	mask(:,2300:end) = NaN;

5872298f   Paul Nicot   maj pour fonction...
40
elseif cam == 'guysB'
a48ac2ca   Paul Nicot   mfile pour le tra...
41
42
	lon1 = -76.8429; lon2 = -75.7753;
	lat1 = 72.5978; lat2 = 72.8511;
79657b4e   Paul Nicot   prise en compte d...
43
	mask = ones(5517,3979);
a48ac2ca   Paul Nicot   mfile pour le tra...
44
45
	mask(:,1:1500) = NaN;
	mask(:,2300:end) = NaN;
79657b4e   Paul Nicot   prise en compte d...
46
    imfiles = dir([rootdir,cam,'/original/2016*_rotated.jpg']);
a48ac2ca   Paul Nicot   mfile pour le tra...
47
48
end

79657b4e   Paul Nicot   prise en compte d...
49
%% x/y resolution
a48ac2ca   Paul Nicot   mfile pour le tra...
50
51
52
53
54
55
56
57
58
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);

79657b4e   Paul Nicot   prise en compte d...
59
60

%% data loading
5872298f   Paul Nicot   maj pour fonction...
61
62
load([ rootdir cam '/work/g_rect_' cam '_20160630.mat']);

a48ac2ca   Paul Nicot   mfile pour le tra...
63
64
65
LON = mask.*LON;
LAT = mask.*LAT;

5872298f   Paul Nicot   maj pour fonction...
66
67
x1 = LON(:); x2 = LON(:);
y1 = LAT(:); 
a48ac2ca   Paul Nicot   mfile pour le tra...
68

79657b4e   Paul Nicot   prise en compte d...
69
70
71
72
73
74
75
76
77
78
79
80
81
82
%% 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)
a48ac2ca   Paul Nicot   mfile pour le tra...
83
84
85
    imfilename = imfiles(n).name;
    imNumber   = imfilename(1:15);

5872298f   Paul Nicot   maj pour fonction...
86
    I0 = imread([ rootdir cam '/original/' imfilename ]);
a48ac2ca   Paul Nicot   mfile pour le tra...
87
    I0 = double(I0);
5872298f   Paul Nicot   maj pour fonction...
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
    
    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'];
a48ac2ca   Paul Nicot   mfile pour le tra...
106
107
108
        
	imwrite(I,jpgname,'jpg');
        
5872298f   Paul Nicot   maj pour fonction...
109
110
111
112
113
114
115
116
    % create geotiff file
    R = georasterref();
    R.RasterSize = [size(I)];
    R.Latlim = [lat1 lat2];
    R.Lonlim = [lon1 lon2];
    R.ColumnsStartFrom = 'north';
    
    geotiffwrite(tifname, I, R);
a48ac2ca   Paul Nicot   mfile pour le tra...
117
118
119
120
121
%    end
    
        
    disp([' ',imNumber]);
end