g_rect.m 13.1 KB
Newer Older
Daniel Bourgault's avatar
Daniel Bourgault committed
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
function g_rect
%G_RECT - Main function for georectifying oblique images.
%
% Syntax:  Simply type g_rect at the command line and follow the
%          instructions.
%
% Inputs:
%    The function G_RECT reads an input parameter file that contains all 
%     information required to perfrorm a georectification of a given image. 
%     The format of this parameter file is detailed on-line on the G_RECT 
%     Wiki page.
%
% Outputs:
%    The function G_RECT creates an output file that contains the following 
%    variables: 
%
%        imgFname:      The reference image for the georectification.
%
%        firstImgFname: The first image of a sequence of images to which the
%                       georectification could be applied to. This is really
%                       just a comment. 
%
%        lastImgFname: The last image of a sequence of images to which the
%                      georectification could be applied to. This is really
%                      just a comment.
%
%        LON:          The main matrix that contain the longitude of each 
%                      pixel of the referecne image (imgFname). Note that
%                      if the package is used to rectify lab images this 
%                      matrix rather contains the distance in meters from
%                      a predefined x-origin. 
%                       
%        LAT:          Same as LON but for the latitude.
% 
%        LON0:         A scalar for the longitude of the camera or, in the 
%                      case of a lab setup its cartesian coordinate in meters.  
%                       
%        LAT0:         Same as LON0 but for the latitude.
%            
%        lon_gcp:      A vector containing the longitude of each ground 
%                      control points (GCP). For the lab case this is the 
%                      cartesian coordinate in meters. 
%                       
%        lat_gcp:      Same as lon_gcp for latitude.
%
%        i_gcp:        The horizontal index of the image ground control
%                      points.
%         
%        j_gcp:        The vertical index of the image ground control
%                      points.
%         
%        hfov:         The camera horizontal field of view [degree].
%
%        phi:          Camera tilt angle [degree].
%
%        H:            The camera altitude relative to the water [m].
%
%        theta:        View angle clockwise from North [degree].
%
%
%        errGeoFit:    The rms error of the georectified image after
%                      geometrical transformation [m].
%
%        errPolyFit:   The rms error of the georectified image after
%                      geometrical transformation and the polynomial 
%                      correction [m].
%
%        precision:    Calculation can be done in single or double
%                      precisions as defined by the user in the parameter
%                      file. With today's computers this is now obsolete 
%                      and calculations can always be done in double 
%                      precision.
%
%
% Other m-files required: Works best with the m_map package for
%                         visualization.
% Subfunctions: all functions contained within the G_RECT folder.
% MAT-files required: none
% 
% Author: Daniel Bourgault
%         Institut des sciences de la mer de Rimouski
% email: daniel_bourgault@uqar.ca 
% Website: http://demeter.uqar.ca/g_rect/
% February 2013
%

%
% The minimization is repeated nMinimize times where each time a random 
% combination of the initial guesses is chosen within the given
% uncertainties provided by the user. This is becasue the algorithm often 
% converges toward a local minimum. The repetition is used to increase chances
% that the minimum found is a true minimum within the uncertainties provided.  
nMinimize = 50;


%% Read the parameter file
% Count the number of header lines before the ground control points (GCPs)
% The GCPs start right after the variable gcpData is set to true.
display('  ');  
display('  Welcome to g_rect: a package for georectifying oblique images on a flat ocean');  
display('  Authors: Daniel Bourgault and Rich Pawlowicz');  
display('  ');  
inputFname = input('  Enter the filename for the input parameters: ','s');  
fid = fopen(inputFname);
 
nHeaderLine = 0;
gcpData = false;

% Read and execute each line of the parameter file until gcpData = true
% after which the GCP data appear and are read below with the importdata
% function.
while gcpData == false
   eval(fgetl(fid));
   nHeaderLine = nHeaderLine + 1;
end
fclose(fid);

%% Import the GCP data (4 column) at the end of the parameter file
gcp = importdata(inputFname,' ',nHeaderLine);
i_gcp   = gcp.data(:,1);
j_gcp   = gcp.data(:,2);
lon_gcp = gcp.data(:,3);
lat_gcp = gcp.data(:,4);
ngcp    = length(i_gcp);

% Get the image size
imgInfo   = imfinfo(imgFname);
imgWidth  = imgInfo.Width;
imgHeight = imgInfo.Height;

if precision == 'single'
  imgWidth  = single(imgWidth);
  imgHeight = single(imgHeight);
end


%% Display information 
fprintf('\n')
fprintf('  INPUT PARAMETERS\n')
fprintf('    Image filename: (imgFname):........... %s\n',imgFname)
fprintf('    First image: (firstImgFname):......... %s\n',firstImgFname)
fprintf('    Last image: (lastImgFname):........... %s\n',lastImgFname)
fprintf('    Output filename: (outputFname):....... %s\n',outputFname);
fprintf('    Image width (imgWidth):............... %i\n',imgWidth)
fprintf('    Image width (imgHeight):.............. %i\n',imgHeight)
fprintf('    Camera longitude (LON0):.............. %f\n',LON0)
fprintf('    Camera latitude (LAT0):............... %f\n',LAT0)
fprintf('    Principal point offset (ic):.......... %f\n',ic)
fprintf('    Principal point offset (jc):.......... %f\n',jc)
fprintf('    Field of view (hfov):................. %f\n',hfov)
fprintf('    Dip angle (lambda):................... %f\n',lambda)
fprintf('    Tilt angle (phi):..................... %f\n',phi)
fprintf('    Camera altitude (H):.................. %f\n',H)
fprintf('    View angle from North (theta):........ %f\n',theta)
fprintf('    Uncertainty in hfov (dhfov):.......... %f\n',dhfov)
fprintf('    Uncertainty in dip angle (dlambda):... %f\n',dlambda)
fprintf('    Uncertainty in tilt angle (dphi):..... %f\n',dphi)
fprintf('    Uncertainty in altitude (dH):......... %f\n',dH)
fprintf('    Uncertainty in view angle (dtheta):... %f\n',dtheta)
fprintf('    Polynomial order (polyOrder):......... %i\n',polyOrder)
fprintf('    Number of GCPs (ngcp):................ %i\n',ngcp)
fprintf('    Precision (precision):................ %s\n',precision)
fprintf('    Field or lab (field=true; lab=false):. %i\n',field)
fprintf('\n')

% Display the image with GCPs;
% image(imread(imgFname));
imagesc(imread(imgFname));
colormap(gray);
170
hold on
Daniel Bourgault's avatar
Daniel Bourgault committed
171
for i = 1:ngcp
172
173
  plot(i_gcp(i),j_gcp(i),'r.');
  text(i_gcp(i),j_gcp(i),num2str(i),'color','r','horizontalalignment','right');
Daniel Bourgault's avatar
Daniel Bourgault committed
174
175
176
end
title('Ground Control Points','color','r');

177
print -dpng -r300 image1.png
Daniel Bourgault's avatar
Daniel Bourgault committed
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
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363

fprintf('\n')
ok = input('Is it ok to proceed with the rectification (y/n): ','s');
if ok ~= 'y'
  return
  %break
end

%%

nUnknown = 0;
if dhfov   > 0.0; nUnknown = nUnknown+1; end
if dlambda > 0.0; nUnknown = nUnknown+1; end
if dphi    > 0.0; nUnknown = nUnknown+1; end
if dH      > 0.0; nUnknown = nUnknown+1; end
if dtheta  > 0.0; nUnknown = nUnknown+1; end

if nUnknown > ngcp
  fprintf('\n')
  fprintf('WARNING: \n');  
  fprintf('         The number of GCPs is < number of unknown parameters.\n');  
  fprintf('         Program stopped.\n');
  %break
  return
end

% Check for consistencies between number of GCPs and order of the polynomial 
% correction
ngcp = length(i_gcp);
if ngcp < 3*polyOrder
  fprintf('\n')
  fprintf('WARNING: \n');  
  fprintf('         The number of GCPs is inconsistent with the order of the polynomial correction.\n');  
  fprintf('         ngcp should be >= 3*polyOrder.\n');  
  fprintf('         Polynomial correction will not be applied.\n');  
  polyCorrection = false;
else
  polyCorrection = true;
end
if polyOrder == 0
  polyCorrection = false;
end

%% This is the main section for the minimization algorithm

if nUnknown > 0
    
  % Options for the fminsearch function. May be needed for some particular
  % problems but in general the default values should work fine.  
  %options=optimset('MaxFunEvals',100000,'MaxIter',100000);
  %options=optimset('MaxFunEvals',100000,'MaxIter',100000,'TolX',1.d-12,'TolFun',1.d-12);
  options = [];
  
    
  % Only feed the minimization algorithm with the GCPs. xp and yp are the
  % image coordinate of these GCPs.
  xp = i_gcp;
  yp = j_gcp;

  % This is the call to the minimization
  bestErrGeoFit = Inf;
  
  % Save inital guesses in new variables. 
  hfovGuess   = hfov;
  lambdaGuess = lambda;
  phiGuess    = phi;
  HGuess      = H;
  thetaGuess  = theta;
  
  for iMinimize = 1:nMinimize
      
      % First guesses for the minimization
      if iMinimize == 1
          hfov0   = hfov; 
          lambda0 = lambda; 
          phi0    = phi; 
          H0      = H;
          theta0  = theta;
      else
          % Select randomly new initial guesses within the specified
          % uncertainties.
          hfov0   = (hfovGuess - dhfov)     + 2*dhfov*rand(1); 
          lambda0 = (lambdaGuess - dlambda) + 2*dlambda*rand(1); 
          phi0    = (phiGuess - dphi)       + 2*dphi*rand(1); 
          H0      = (HGuess - dH)           + 2*dH*rand(1);
          theta0  = (thetaGuess - dtheta)   + 2*dtheta*rand(1);
      end
  
      % Cretae vector cv0 for the initial guesses. 
      i = 0;
      if dhfov > 0.0
          i = i+1;
          cv0(i) = hfov0;
          theOrder(i) = 1;
      end
      if dlambda > 0.0
          i = i + 1;
          cv0(i) = lambda0;
          theOrder(i) = 2;
      end
      if dphi > 0.0
          i = i + 1;
          cv0(i) = phi0;
          theOrder(i) = 3;
      end
      if dH > 0.0
          i = i + 1;
          cv0(i) = H0;
          theOrder(i) = 4;
      end
      if dtheta > 0.0
          i = i + 1;
          cv0(i) = theta0;
          theOrder(i) = 5;
      end

      [cv, errGeoFit] = fminsearch('g_error_geofit',cv0,options, ...
                                  imgWidth,imgHeight,xp,yp,ic,jc,...
                                  hfov,lambda,phi,H,theta,...
                                  hfov0,lambda0,phi0,H0,theta0,...
                                  hfovGuess,lambdaGuess,phiGuess,HGuess,thetaGuess,...
                                  dhfov,dlambda,dphi,dH,dtheta,...
                                  LON0,LAT0,...
                                  i_gcp,j_gcp,lon_gcp,lat_gcp,...
                                  theOrder,field);

                              
      if errGeoFit < bestErrGeoFit
          bestErrGeoFit = errGeoFit;
          cvBest = cv;
      end
      
      fprintf('\n')
      fprintf('  Iteration # (iMinimize):                       %i\n',iMinimize);
      fprintf('  Max. number of iteration (nMinimize):          %i\n',nMinimize);
      fprintf('  RSM error (m)  for this iteration (errGeoFit): %f\n',errGeoFit);
      fprintf('  Best RSM error (m) so far (bestErrGeoFit):     %f\n',bestErrGeoFit);
      
  end
  
  for i = 1:length(theOrder)
    if theOrder(i) == 1; hfov   = cvBest(i); end
    if theOrder(i) == 2; lambda = cvBest(i); end
    if theOrder(i) == 3; phi    = cvBest(i); end
    if theOrder(i) == 4; H      = cvBest(i); end
    if theOrder(i) == 5; theta  = cvBest(i); end
  end
  
  fprintf('\n')
  fprintf('PARAMETERS AFTER GEOMETRICAL RECTIFICATION \n')
  fprintf('  Field of view (hfov):            %f\n',hfov)
  fprintf('  Dip angle (lambda):              %f\n',lambda)
  fprintf('  Tilt angle (phi):                %f\n',phi)
  fprintf('  Camera altitude (H):             %f\n',H)
  fprintf('  View angle from North (theta):   %f\n',theta)
  fprintf('\n')

  if length(theOrder) > 1
    fprintf('The rms error after geometrical correction (m): %f\n',bestErrGeoFit);
  end
  
end

%%  

% Now construct the matrices LON and LAT for the entire image using the 
% camera parameters found by minimization just above.

% Camera coordinate of all pixels
xp = repmat([1:imgWidth]',1,imgHeight);
yp = repmat([1:imgHeight],imgWidth,1);

% Transform camera coordinate to ground coordinate.
[LON LAT] = g_pix2ll(xp,yp,imgWidth,imgHeight,ic,jc,...
                     hfov,lambda,phi,theta,H,LON0,LAT0,field);


%% Apply polynomial correction if requested.
if polyCorrection == true
  [LON LAT errPolyFit] = g_poly(LON,LAT,LON0,LAT0,i_gcp,j_gcp,lon_gcp,lat_gcp,polyOrder,field);
  fprintf('The rms error after polynomial stretching (m):  %f\n',errPolyFit)
else
  errPolyFit = NaN;
end
%%

364
365
366
% Compute the effective resolution
Delta = g_res(LON, LAT, field);

Daniel Bourgault's avatar
Daniel Bourgault committed
367
368
369
370
371
372
373
374
375
376
fprintf('\n')
fprintf('Saving rectification file in: %s\n',outputFname);

save(outputFname,'imgFname','firstImgFname','lastImgFname',...
                 'LON','LAT',...
                 'LON0','LAT0',...
                 'lon_gcp','lat_gcp',...
                 'i_gcp','j_gcp',...
                 'hfov','lambda','phi','H','theta',...
                 'errGeoFit','errPolyFit',...
377
                 'precision','Delta');
Daniel Bourgault's avatar
Daniel Bourgault committed
378
379
380
381
382
383
384
385
386
387
388
389
390
391

clear LON LAT

fprintf('\n')
fprintf('Making figure\n');

figure;  
if field
    g_viz_field(imgFname,outputFname);
else
    g_viz_lab(imgFname,outputFname);
end

print -dpng image2.png