g_rect.m 14.5 KB
Newer Older
Daniel Bourgault's avatar
Daniel Bourgault committed
1
2
3
4
5
6
7
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:
8
9
10
%    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
Daniel Bourgault's avatar
Daniel Bourgault committed
11
12
13
%     Wiki page.
%
% Outputs:
14
15
%    The function G_RECT creates an output file that contains the following
%    variables:
Daniel Bourgault's avatar
Daniel Bourgault committed
16
17
18
19
20
%
%        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
21
%                       just a comment.
Daniel Bourgault's avatar
Daniel Bourgault committed
22
23
24
25
26
%
%        lastImgFname: The last image of a sequence of images to which the
%                      georectification could be applied to. This is really
%                      just a comment.
%
27
%        LON:          The main matrix that contain the longitude of each
Daniel Bourgault's avatar
Daniel Bourgault committed
28
%                      pixel of the referecne image (imgFname). Note that
29
%                      if the package is used to rectify lab images this
Daniel Bourgault's avatar
Daniel Bourgault committed
30
%                      matrix rather contains the distance in meters from
31
32
%                      a predefined x-origin.
%
Daniel Bourgault's avatar
Daniel Bourgault committed
33
%        LAT:          Same as LON but for the latitude.
34
35
36
37
%
%        LON0:         A scalar for the longitude of the camera or, in the
%                      case of a lab setup its cartesian coordinate in meters.
%
Daniel Bourgault's avatar
Daniel Bourgault committed
38
%        LAT0:         Same as LON0 but for the latitude.
39
40
41
42
43
44
45
46
47
48
49
50
51
%
%        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.
%
Daniel Bourgault's avatar
Daniel Bourgault committed
52
53
%        lat_gcp:      Same as lon_gcp for latitude.
%
54
55
56
%        h_gcp:        The elevation in meter of the GCPs. The elevation is
%                      0 m if taken at water level.
%
Daniel Bourgault's avatar
Daniel Bourgault committed
57
58
%        i_gcp:        The horizontal index of the image ground control
%                      points.
59
%
Daniel Bourgault's avatar
Daniel Bourgault committed
60
61
%        j_gcp:        The vertical index of the image ground control
%                      points.
62
%
Daniel Bourgault's avatar
Daniel Bourgault committed
63
64
65
66
%        hfov:         The camera horizontal field of view [degree].
%
%        phi:          Camera tilt angle [degree].
%
67
68
%        lambda:       Camera dip angle [degree].
%
Daniel Bourgault's avatar
Daniel Bourgault committed
69
70
71
72
73
74
75
76
%        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
77
%                      geometrical transformation and the polynomial
Daniel Bourgault's avatar
Daniel Bourgault committed
78
79
80
81
%                      correction [m].
%
%        precision:    Calculation can be done in single or double
%                      precisions as defined by the user in the parameter
82
83
%                      file. With today's computers this is now obsolete
%                      and calculations can always be done in double
Daniel Bourgault's avatar
Daniel Bourgault committed
84
85
%                      precision.
%
86
%         Delta:       The effective resolution of the georectified image.
Daniel Bourgault's avatar
Daniel Bourgault committed
87
88
89
90
%
% Other m-files required: Works best with the m_map package for
%                         visualization.
% Subfunctions: all functions contained within the G_RECT folder.
91
%
Daniel Bourgault's avatar
Daniel Bourgault committed
92
93
% Author: Daniel Bourgault
%         Institut des sciences de la mer de Rimouski
94
% email: daniel_bourgault@uqar.ca
95
% Website: https://srwebpolr01.uqar.ca/g_rect/
Daniel Bourgault's avatar
Daniel Bourgault committed
96
97
98
% February 2013
%

99
100
%%
% The minimization is repeated nMinimize times where each time a random
Daniel Bourgault's avatar
Daniel Bourgault committed
101
% combination of the initial guesses is chosen within the given
102
% uncertainties provided by the user. This is becasue the algorithm often
Daniel Bourgault's avatar
Daniel Bourgault committed
103
% converges toward a local minimum. The repetition is used to increase chances
104
% that the minimum found is a true minimum within the uncertainties provided.
Daniel Bourgault's avatar
Daniel Bourgault committed
105
106
107
108
109
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.
110
111
112
113
114
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');
Daniel Bourgault's avatar
Daniel Bourgault committed
115
fid = fopen(inputFname);
116
117
118

nHeaderLine  = 0;
gcpData      = false;
Daniel Bourgault's avatar
Daniel Bourgault committed
119
120
121
122
123

% 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
124
125
    eval(fgetl(fid));
    nHeaderLine = nHeaderLine + 1;
Daniel Bourgault's avatar
Daniel Bourgault committed
126
127
128
end
fclose(fid);

129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
%% 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);
h_gcp    = gcp.data(:,5);

%%
% 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
ngcp = length(i_gcp);
Daniel Bourgault's avatar
Daniel Bourgault committed
152
153
154
155
156
157
158

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

if precision == 'single'
159
160
    imgWidth  = single(imgWidth);
    imgHeight = single(imgHeight);
Daniel Bourgault's avatar
Daniel Bourgault committed
161
162
end

163
%% Display information
Daniel Bourgault's avatar
Daniel Bourgault committed
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
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);
196
hold on
Daniel Bourgault's avatar
Daniel Bourgault committed
197
for i = 1:ngcp
198
199
200
201
202
    plot(i_gcp(i),j_gcp(i),'r.');
    text(i_gcp(i),j_gcp(i),[num2str(i),'(',num2str(h_gcp(i)),')'],...
        'color','r',...
        'horizontalalignment','right',...
        'fontsize',16);
Daniel Bourgault's avatar
Daniel Bourgault committed
203
end
204
title('Ground Control Points Number (elevation in meters)','color','r');
Daniel Bourgault's avatar
Daniel Bourgault committed
205

206
print -dpng -r300 image1.png
Daniel Bourgault's avatar
Daniel Bourgault committed
207
208
209
210

fprintf('\n')
ok = input('Is it ok to proceed with the rectification (y/n): ','s');
if ok ~= 'y'
211
    return
Daniel Bourgault's avatar
Daniel Bourgault committed
212
213
214
215
216
217
218
219
220
221
222
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
223
224
225
226
227
    fprintf('\n')
    fprintf('WARNING: \n');
    fprintf('         The number of GCPs is < number of unknown parameters.\n');
    fprintf('         Program stopped.\n');
    return
Daniel Bourgault's avatar
Daniel Bourgault committed
228
229
end

230
% Check for consistencies between number of GCPs and order of the polynomial
Daniel Bourgault's avatar
Daniel Bourgault committed
231
232
233
% correction
ngcp = length(i_gcp);
if ngcp < 3*polyOrder
234
235
236
237
238
239
    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;
Daniel Bourgault's avatar
Daniel Bourgault committed
240
else
241
    polyCorrection = true;
Daniel Bourgault's avatar
Daniel Bourgault committed
242
243
end
if polyOrder == 0
244
    polyCorrection = false;
Daniel Bourgault's avatar
Daniel Bourgault committed
245
246
247
248
249
250
end

%% This is the main section for the minimization algorithm

if nUnknown > 0
    
251
252
253
254
255
256
257
258
259
260
    % 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;
Daniel Bourgault's avatar
Daniel Bourgault committed
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
    % 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_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
Daniel Bourgault's avatar
Daniel Bourgault committed
362
363
end

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

%%
Daniel Bourgault's avatar
Daniel Bourgault committed
368

369
% Now construct the matrices LON and LAT for the entire image using the
Daniel Bourgault's avatar
Daniel Bourgault committed
370
371
372
373
374
375
376
377
378
379
380
381
382
% 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
383
384
    [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)
Daniel Bourgault's avatar
Daniel Bourgault committed
385
else
386
    errPolyFit = NaN;
Daniel Bourgault's avatar
Daniel Bourgault committed
387
388
389
end
%%

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

Daniel Bourgault's avatar
Daniel Bourgault committed
393
394
395
396
397
398
fprintf('\n')
fprintf('Saving rectification file in: %s\n',outputFname);

save(outputFname,'imgFname','firstImgFname','lastImgFname',...
                 'LON','LAT',...
                 'LON0','LAT0',...
399
400
                 'lon_gcp0','lat_gcp0',...
                 'lon_gcp','lat_gcp','h_gcp',...
Daniel Bourgault's avatar
Daniel Bourgault committed
401
402
403
                 'i_gcp','j_gcp',...
                 'hfov','lambda','phi','H','theta',...
                 'errGeoFit','errPolyFit',...
404
                 'precision','Delta');
Daniel Bourgault's avatar
Daniel Bourgault committed
405
406
407
408
409
410

clear LON LAT

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

411
figure;
Daniel Bourgault's avatar
Daniel Bourgault committed
412
413
414
415
416
417
if field
    g_viz_field(imgFname,outputFname);
else
    g_viz_lab(imgFname,outputFname);
end

418
print -dpng image2.png