-
Notifications
You must be signed in to change notification settings - Fork 22
/
Copy pathdemoengabreen.m
222 lines (165 loc) · 7.17 KB
/
demoengabreen.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
%% Feature tracking example
%
% This is a complete example of feature tracking on Engabreen.
%
% # Load images & data
% # Use GCPs to determine camera view direction and lens distortion
% parameters of image
% # track stable rock features to determine camera shake and infer view direction
% of image B
% # Pre-process DEM by filling crevasses.
% # Track ice motion between images
% # Georeference tracked points and calculate real world velocities.
%
%
close all
%% Setup file locations and load images & data
%
%
idA = 8902; idB = 8937; % image ids (/file numbers)
datafolder = 'demos';
fA = fullfile(datafolder,sprintf('IMG_%4.0f.jpg',idA));
fB = fullfile(datafolder,sprintf('IMG_%4.0f.jpg',idB));
%load images:
A = imread(fA);
B = imread(fB);
metaA = imfinfo(fA);tA = datenum(metaA.DateTime,'yyyy:mm:dd HH:MM:SS');
metaB = imfinfo(fB);tB = datenum(metaB.DateTime,'yyyy:mm:dd HH:MM:SS');
dem = load(fullfile(datafolder,'dem')); %load DEM
gcpA = load(fullfile(datafolder,'gcp8902.txt'));%load ground control points for image A
%% Determine camera parameters for image A
%
% # Initial crude guess at camera parameters
% # Use GCPs to optimize camera parameters
%
%calculate focal length in pixel units:
FocalLength = 30; %mm (can also be found here: metaA.DigitalCamera.FocalLength)
SensorSize = [22.0 14.7]; %mm: http://www.cnet.com/products/canon-eos-rebel-t3/specs/
imgsz = size(A);
f = imgsz([2 1]).*(FocalLength./SensorSize);
%known camera location:
cameralocation = [446722.0 7396671.0 770.0];
%crude estimate of look direction.
camA = camera(cameralocation,size(A),[200 0 0]*pi/180,f); %loooking west
%Use GCPs to optimize the following camera parameters:
%view dir, focal lengths, and a simple radial distortion model
[camA,rmse,aic] = camA.optimizecam(gcpA(:,1:3),gcpA(:,4:5),'00000111110010000000');
fprintf('reprojectionerror = %3.1fpx AIC:%4.0f\n',rmse,aic)
%Visually compare the projection of the GCPs with the pixel coords:
figure
axes('position',[0 .1 1 .8]); hold on
image(A)
axis equal off ij tight
hold on
uv = camA.project(gcpA(:,1:3));
h = plot(gcpA(:,4),gcpA(:,5),'+',uv(:,1),uv(:,2),'rx');
legend(h,'UV of GCP','projection of GCPs','location','southoutside')
title(sprintf('Projection of ground control points. RMSE = %.1fpx',rmse))
%% Determine view direction of camera B.
%
% # find movement of rock features between images A and B
% # determine camera B by pertubing viewdir of camera A.
% First get an approximate estimate of the image shift using a single large
% template
[duoffset,dvoffset] = templatematch(A,B,3000,995,'templatewidth',261,'searchwidth',400,'supersample',0.5,'showprogress',false)
% Get a whole bunch of image shift estimates using a grid of probe points.
% Having multiple shift estimates will allow us to determine camera
% rotation.
[pu,pv] = meshgrid(200:700:4000,100:400:1000);
pu = pu(:); pv = pv(:)+pu/10;
[du,dv,C] = templatematch(A,B,pu,pv,'templatewidth',61,'searchwidth',81,'supersample',3,'initialdu',duoffset,'initialdv',dvoffset);
% Determine camera rotation between A and B from the set of image
% shifts.
% find 3d coords consistent with the 2d pixel coords in points.
xyz = camA.invproject([pu pv]);
% the projection of xyz has to match the shifted coords in points+dxy:
[camB,rmse] = camA.optimizecam(xyz,[pu+du pv+dv],'00000111000000000000'); %optimize 3 view direction angles to determine camera B.
rmse
%quantify the shift between A and B in terms of an delta angle.
DeltaViewDirection = (camB.viewdir-camA.viewdir)*180/pi
%% Generate a set of points to be tracked between images
%
% # Generate a regular grid of candidate points in world coordinates.
% # Cull the set of candidate points to those that are visible and glaciated
%
%
% The viewshed is all the points of the dem that are visible from the
% camera location. They may not be in the field of view of the lens.
dem.visible = voxelviewshed(dem.X,dem.Y,dem.filled,camA.xyz);
%Make a regular 50 m grid of points we would like to track from image A
[XA,YA] = meshgrid(min(dem.x):50:max(dem.x),min(dem.y):50:max(dem.y));
ZA = interp2(dem.X,dem.Y,dem.filled,XA,YA);
%Figure out which pixel coordinates they correspond to:
[uvA,~,inframe] = camA.project([XA(:) YA(:) ZA(:)]); %where would the candidate points be in image A
%Insert nans where we do not want to track:
keepers = double(dem.visible&dem.mask); %visible & glaciated dem points
keepers = filter2(ones(11)/(11^2),keepers); %throw away points close to the edge of visibility
keepers = interp2(dem.X,dem.Y,keepers,XA(:),YA(:))>.99; %which candidate points fullfill the criteria.
uvA(~(keepers&inframe)) = nan;
%% Track points between images.
% calculate where points would be in image B if no ice motion
% ( i.e. accounting only for camera shake)
camshake = camB.project(camA.invproject(uvA))-uvA;
options = [];
options.pu = uvA(:,1);
options.pv = uvA(:,2);
options.method = 'OC';
options.showprogress = true;
options.searchwidth = 81;
options.templatewidth = 21;
options.supersample = 2; %supersample the input images for better subpixel estimation
options.initialdu = camshake(:,1);
options.initialdv = camshake(:,2);
[du,dv,C,Cnoise] = templatematch(A,B,options);
uvB = uvA+[du dv];
signal2noise = C./Cnoise;
%% Georeference tracked points
% ... and calculate velocities
xyzA = camA.invproject(uvA,dem.X,dem.Y,dem.filled); % has to be recalculated because uvA has been rounded.
xyzB = camB.invproject(uvB,dem.X,dem.Y,dem.filled-dem.mask*22.75*(tB-tA)/365); % impose a thinning of the DEM of 23m/yr between images.
V = (xyzB-xyzA)./(tB-tA); % 3d velocity.
Vx = reshape(V(:,1),size(XA));
Vy = reshape(V(:,2),size(YA));
Vz = reshape(V(:,3),size(ZA));
figure;
showimg(dem.x,dem.y,dem.rgb);
hold on
Vn = sqrt(sum(V(:,1:2).^2,2));
keep = signal2noise>2 & C>.7;
alphawarp(XA,YA,sqrt(Vx.^2+Vy.^2))
quiver(xyzA(keep,1),xyzA(keep,2),V(keep,1)./Vn(keep),V(keep,2)./Vn(keep),.2,'k')
caxis([0 1])
colormap hot
hcb = colorbar('southoutside');
plot(camA.xyz(1),camA.xyz(2),'r+')
title('Velocity in metres per day')
%% Project velocity onto downhill slope direction
% ----
% The largest error in the velocities will along the view direction vector.
% By projecting to the slope direction we strongly suppress errors arising
% from this.
[gradX,gradY] = gradient(dem.filled,dem.X(2,2)-dem.X(1,1),dem.Y(2,2)-dem.Y(1,1));
gradN = sqrt(gradX.^2+gradY.^2);
gradX = -gradX./gradN;gradY = -gradY./gradN;
gradX = interp2(dem.X,dem.Y,gradX,XA,YA);
gradY = interp2(dem.X,dem.Y,gradY,XA,YA);
Vgn = Vx.*gradX+Vy.*gradY;%Velocity along glacier
Vacross = Vx.*gradY-Vy.*gradX; %Velocity across glacier
keep = reshape(keep,size(XA));
%We do not trust regions with large across glacier flow. We may get large errors where the motion is
%in and out of the frame. I.e. in places where the glacier surface is
%viewed very obliquely. In those places we do not have a sufficiently good
%view to calculate velocities. We apply a filter to remove these.
keep = keep&(abs(Vgn)>abs(Vacross));
close all
figure
showimg(dem.x,dem.y,dem.rgb);
axis equal xy off tight
hold on
alphawarp(XA,YA,Vgn,keep*.7)
quiver(XA(keep),YA(keep),gradX(keep),gradY(keep),.2,'k')
caxis([0 1])
colormap hot
hcb = colorbar('southoutside');
plot(camA.xyz(1),camA.xyz(2),'r+')
title('Velocity along slope direction in metres per day')