g_rect.m 15.2 KB
Newer Older
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
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
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
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.
%
%        lon0_gcp:     A vector containing the longitude of each ground
%                      control points (GCP). For the lab case this is the
%                      cartesian coordinate in meters. These GCPs may have
%                      elevations.
%
%        lat0_gcp:     Same as lon_gcp for latitude.
%
%        lon_gcp:      A vector containing the longitude of each ground
%                      control points (GCP) projected onto the water level
%                      i.e. at 0 m of elevation. For the lab case this is 
%                      the cartesian coordinate in meters.
%
%        lat_gcp:      Same as lon_gcp for latitude.
%
%        h_gcp:        The elevation in meter of the GCPs. The elevation is
%                      0 m if taken at water level.
%
%        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].
%
%        lambda:       Camera dip 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.
%
%         Delta:       The effective resolution of the georectified image.
%
% Other m-files required: Works best with the m_map package for
%                         visualization.
% Subfunctions: all functions contained within the G_RECT folder.
%
% Author: Daniel Bourgault
%         Institut des sciences de la mer de Rimouski
% email: daniel_bourgault@uqar.ca
% Website: https://srwebpolr01.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 at the end of the parameter file
gcp      = importdata(inputFname,' ',nHeaderLine);
i_gcp    = gcp.data(:,1);
j_gcp    = gcp.data(:,2);
lon_gcp0 = gcp.data(:,3);
lat_gcp0 = gcp.data(:,4);

[ngcp n_column] = size(gcp.data)

% If there are 5 columns it means that elevation are provided
if n_column == 5
    h_gcp    = gcp.data(:,5);
else % otherwise elevation are considered to be zero
    h_gcp(1:ngcp) = 0.0;
end

%%
% Check if the elevation of the GCPs are not too high and above
% a certain fraction (gamma) of the camera height. If so, stop.
gamma = 0.75;
i_bad = find(h_gcp > gamma*(H+dH));
if length(i_bad) > 0
    display([' ']);
    display(['  WARNING:']);
    for i = 1:length(i_bad)
        display(['      The elevation of GCP #',num2str(i_bad(i)),' is greater than ',num2str(gamma),'*(H+dH).']);
    end
    display(['  FIX AND RERUN.']);
    return
end

% 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);
hold on
for i = 1:ngcp
    plot(i_gcp(i),j_gcp(i),'r.');
    text(i_gcp(i),j_gcp(i),[' ',num2str(i),'(',num2str(h_gcp(i)),')'],...
        'color','r',...
        'horizontalalignment','left',...
        'fontsize',10);
end
title('Ground Control Points Number (elevation in meters)','color','r');
daspect([1 1 1]);
print('-dpng','-r300',[imgFname(1:end-4),'_GCP.png']);

fprintf('\n')
ok = input('Is it ok to proceed with the rectification (y/n): ','s');
if ok ~= 'y'
    return
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');
    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
        
        % Create 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_gcp0,lat_gcp0,h_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

% Project the GCP that have elevation.
[lon_gcp,lat_gcp] = g_proj_GCP(LON0,LAT0,H,lon_gcp0,lat_gcp0,h_gcp,field);

%%

% 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
%%

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

fprintf('\n')
fprintf('Saving rectification file in: %s\n',outputFname);

save(outputFname,'imgFname','firstImgFname','lastImgFname',...
                 'LON','LAT',...
                 'LON0','LAT0',...
                 'lon_gcp0','lat_gcp0',...
                 'lon_gcp','lat_gcp','h_gcp',...
                 'i_gcp','j_gcp',...
                 'hfov','lambda','phi','H','theta',...
                 'errGeoFit','errPolyFit',...
                 'precision','Delta');

clear LON LAT

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

if field
    g_viz_field(imgFname,outputFname);
else
    g_viz_lab(imgFname,outputFname);
end
print('-dpng','-r300',[imgFname(1:end-4),'_grect.png']);