Commit 1eda8cb8dd5702384236e4b2efd976d7dd49a927

Authored by Daniel Bourgault
1 parent 729a68ae
Exists in master

Change the nameing convention of the figures produced (_GCP and _grect).

Showing 1 changed file with 423 additions and 424 deletions   Show diff stats
src/g_rect/g_rect.m 100755 → 100644
1   -function g_rect
2   -%G_RECT - Main function for georectifying oblique images.
3   -%
4   -% Syntax: Simply type g_rect at the command line and follow the
5   -% instructions.
6   -%
7   -% Inputs:
8   -% The function G_RECT reads an input parameter file that contains all
9   -% information required to perfrorm a georectification of a given image.
10   -% The format of this parameter file is detailed on-line on the G_RECT
11   -% Wiki page.
12   -%
13   -% Outputs:
14   -% The function G_RECT creates an output file that contains the following
15   -% variables:
16   -%
17   -% imgFname: The reference image for the georectification.
18   -%
19   -% firstImgFname: The first image of a sequence of images to which the
20   -% georectification could be applied to. This is really
21   -% just a comment.
22   -%
23   -% lastImgFname: The last image of a sequence of images to which the
24   -% georectification could be applied to. This is really
25   -% just a comment.
26   -%
27   -% LON: The main matrix that contain the longitude of each
28   -% pixel of the referecne image (imgFname). Note that
29   -% if the package is used to rectify lab images this
30   -% matrix rather contains the distance in meters from
31   -% a predefined x-origin.
32   -%
33   -% LAT: Same as LON but for the latitude.
34   -%
35   -% LON0: A scalar for the longitude of the camera or, in the
36   -% case of a lab setup its cartesian coordinate in meters.
37   -%
38   -% LAT0: Same as LON0 but for the latitude.
39   -%
40   -% lon0_gcp: A vector containing the longitude of each ground
41   -% control points (GCP). For the lab case this is the
42   -% cartesian coordinate in meters. These GCPs may have
43   -% elevations.
44   -%
45   -% lat0_gcp: Same as lon_gcp for latitude.
46   -%
47   -% lon_gcp: A vector containing the longitude of each ground
48   -% control points (GCP) projected onto the water level
49   -% i.e. at 0 m of elevation. For the lab case this is
50   -% the cartesian coordinate in meters.
51   -%
52   -% lat_gcp: Same as lon_gcp for latitude.
53   -%
54   -% h_gcp: The elevation in meter of the GCPs. The elevation is
55   -% 0 m if taken at water level.
56   -%
57   -% i_gcp: The horizontal index of the image ground control
58   -% points.
59   -%
60   -% j_gcp: The vertical index of the image ground control
61   -% points.
62   -%
63   -% hfov: The camera horizontal field of view [degree].
64   -%
65   -% phi: Camera tilt angle [degree].
66   -%
67   -% lambda: Camera dip angle [degree].
68   -%
69   -% H: The camera altitude relative to the water [m].
70   -%
71   -% theta: View angle clockwise from North [degree].
72   -%
73   -% errGeoFit: The rms error of the georectified image after
74   -% geometrical transformation [m].
75   -%
76   -% errPolyFit: The rms error of the georectified image after
77   -% geometrical transformation and the polynomial
78   -% correction [m].
79   -%
80   -% precision: Calculation can be done in single or double
81   -% precisions as defined by the user in the parameter
82   -% file. With today's computers this is now obsolete
83   -% and calculations can always be done in double
84   -% precision.
85   -%
86   -% Delta: The effective resolution of the georectified image.
87   -%
88   -% Other m-files required: Works best with the m_map package for
89   -% visualization.
90   -% Subfunctions: all functions contained within the G_RECT folder.
91   -%
92   -% Author: Daniel Bourgault
93   -% Institut des sciences de la mer de Rimouski
94   -% email: daniel_bourgault@uqar.ca
95   -% Website: https://srwebpolr01.uqar.ca/g_rect/
96   -% February 2013
97   -%
98   -
99   -%%
100   -% The minimization is repeated nMinimize times where each time a random
101   -% combination of the initial guesses is chosen within the given
102   -% uncertainties provided by the user. This is becasue the algorithm often
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.
105   -nMinimize = 50;
106   -
107   -%% Read the parameter file
108   -% Count the number of header lines before the ground control points (GCPs)
109   -% The GCPs start right after the variable gcpData is set to true.
110   -display(' ');
111   -display(' Welcome to g_rect: a package for georectifying oblique images on a flat ocean');
112   -display(' Authors: Daniel Bourgault and Rich Pawlowicz');
113   -display(' ');
114   -inputFname = input(' Enter the filename for the input parameters: ','s');
115   -fid = fopen(inputFname);
116   -
117   -nHeaderLine = 0;
118   -gcpData = false;
119   -
120   -% Read and execute each line of the parameter file until gcpData = true
121   -% after which the GCP data appear and are read below with the importdata
122   -% function.
123   -while gcpData == false
124   - eval(fgetl(fid));
125   - nHeaderLine = nHeaderLine + 1;
126   -end
127   -fclose(fid);
128   -
129   -%% Import the GCP data at the end of the parameter file
130   -gcp = importdata(inputFname,' ',nHeaderLine);
131   -i_gcp = gcp.data(:,1);
132   -j_gcp = gcp.data(:,2);
133   -lon_gcp0 = gcp.data(:,3);
134   -lat_gcp0 = gcp.data(:,4);
135   -
136   -[ngcp n_column] = size(gcp.data)
137   -
138   -% If there are 5 columns it means that elevation are provided
139   -if n_column == 5
140   - h_gcp = gcp.data(:,5);
141   -else % otherwise elevation are considered to be zero
142   - h_gcp(1:ngcp) = 0.0;
143   -end
144   -
145   -%%
146   -% Check if the elevation of the GCPs are not too high and above
147   -% a certain fraction (gamma) of the camera height. If so, stop.
148   -gamma = 0.75;
149   -i_bad = find(h_gcp > gamma*(H+dH));
150   -if length(i_bad) > 0
151   - display([' ']);
152   - display([' WARNING:']);
153   - for i = 1:length(i_bad)
154   - display([' The elevation of GCP #',num2str(i_bad(i)),' is greater than ',num2str(gamma),'*(H+dH).']);
155   - end
156   - display([' FIX AND RERUN.']);
157   - return
158   -end
159   -
160   -% Get the image size
161   -imgInfo = imfinfo(imgFname);
162   -imgWidth = imgInfo.Width;
163   -imgHeight = imgInfo.Height;
164   -
165   -if precision == 'single'
166   - imgWidth = single(imgWidth);
167   - imgHeight = single(imgHeight);
168   -end
169   -
170   -%% Display information
171   -fprintf('\n')
172   -fprintf(' INPUT PARAMETERS\n')
173   -fprintf(' Image filename: (imgFname):........... %s\n',imgFname)
174   -fprintf(' First image: (firstImgFname):......... %s\n',firstImgFname)
175   -fprintf(' Last image: (lastImgFname):........... %s\n',lastImgFname)
176   -fprintf(' Output filename: (outputFname):....... %s\n',outputFname);
177   -fprintf(' Image width (imgWidth):............... %i\n',imgWidth)
178   -fprintf(' Image width (imgHeight):.............. %i\n',imgHeight)
179   -fprintf(' Camera longitude (LON0):.............. %f\n',LON0)
180   -fprintf(' Camera latitude (LAT0):............... %f\n',LAT0)
181   -fprintf(' Principal point offset (ic):.......... %f\n',ic)
182   -fprintf(' Principal point offset (jc):.......... %f\n',jc)
183   -fprintf(' Field of view (hfov):................. %f\n',hfov)
184   -fprintf(' Dip angle (lambda):................... %f\n',lambda)
185   -fprintf(' Tilt angle (phi):..................... %f\n',phi)
186   -fprintf(' Camera altitude (H):.................. %f\n',H)
187   -fprintf(' View angle from North (theta):........ %f\n',theta)
188   -fprintf(' Uncertainty in hfov (dhfov):.......... %f\n',dhfov)
189   -fprintf(' Uncertainty in dip angle (dlambda):... %f\n',dlambda)
190   -fprintf(' Uncertainty in tilt angle (dphi):..... %f\n',dphi)
191   -fprintf(' Uncertainty in altitude (dH):......... %f\n',dH)
192   -fprintf(' Uncertainty in view angle (dtheta):... %f\n',dtheta)
193   -fprintf(' Polynomial order (polyOrder):......... %i\n',polyOrder)
194   -fprintf(' Number of GCPs (ngcp):................ %i\n',ngcp)
195   -fprintf(' Precision (precision):................ %s\n',precision)
196   -fprintf(' Field or lab (field=true; lab=false):. %i\n',field)
197   -fprintf('\n')
198   -
199   -% Display the image with GCPs;
200   -% image(imread(imgFname));
201   -imagesc(imread(imgFname));
202   -colormap(gray);
203   -hold on
204   -for i = 1:ngcp
205   - plot(i_gcp(i),j_gcp(i),'r.');
206   - text(i_gcp(i),j_gcp(i),[num2str(i),'(',num2str(h_gcp(i)),')'],...
207   - 'color','r',...
208   - 'horizontalalignment','right',...
209   - 'fontsize',16);
210   -end
211   -title('Ground Control Points Number (elevation in meters)','color','r');
212   -
213   -print -dpng -r300 image1.png
214   -
215   -fprintf('\n')
216   -ok = input('Is it ok to proceed with the rectification (y/n): ','s');
217   -if ok ~= 'y'
218   - return
219   -end
220   -
221   -%%
222   -nUnknown = 0;
223   -if dhfov > 0.0; nUnknown = nUnknown+1; end
224   -if dlambda > 0.0; nUnknown = nUnknown+1; end
225   -if dphi > 0.0; nUnknown = nUnknown+1; end
226   -if dH > 0.0; nUnknown = nUnknown+1; end
227   -if dtheta > 0.0; nUnknown = nUnknown+1; end
228   -
229   -if nUnknown > ngcp
230   - fprintf('\n')
231   - fprintf('WARNING: \n');
232   - fprintf(' The number of GCPs is < number of unknown parameters.\n');
233   - fprintf(' Program stopped.\n');
234   - return
235   -end
236   -
237   -% Check for consistencies between number of GCPs and order of the polynomial
238   -% correction
239   -ngcp = length(i_gcp);
240   -if ngcp < 3*polyOrder
241   - fprintf('\n')
242   - fprintf('WARNING: \n');
243   - fprintf(' The number of GCPs is inconsistent with the order of the polynomial correction.\n');
244   - fprintf(' ngcp should be >= 3*polyOrder.\n');
245   - fprintf(' Polynomial correction will not be applied.\n');
246   - polyCorrection = false;
247   -else
248   - polyCorrection = true;
249   -end
250   -if polyOrder == 0
251   - polyCorrection = false;
252   -end
253   -
254   -%% This is the main section for the minimization algorithm
255   -
256   -if nUnknown > 0
257   -
258   - % Options for the fminsearch function. May be needed for some particular
259   - % problems but in general the default values should work fine.
260   - %options=optimset('MaxFunEvals',100000,'MaxIter',100000);
261   - %options=optimset('MaxFunEvals',100000,'MaxIter',100000,'TolX',1.d-12,'TolFun',1.d-12);
262   - options = [];
263   -
264   - % Only feed the minimization algorithm with the GCPs. xp and yp are the
265   - % image coordinate of these GCPs.
266   - xp = i_gcp;
267   - yp = j_gcp;
268   -
269   - % This is the call to the minimization
270   - bestErrGeoFit = Inf;
271   -
272   - % Save inital guesses in new variables.
273   - hfovGuess = hfov;
274   - lambdaGuess = lambda;
275   - phiGuess = phi;
276   - HGuess = H;
277   - thetaGuess = theta;
278   -
279   - for iMinimize = 1:nMinimize
280   -
281   - % First guesses for the minimization
282   - if iMinimize == 1
283   - hfov0 = hfov;
284   - lambda0 = lambda;
285   - phi0 = phi;
286   - H0 = H;
287   - theta0 = theta;
288   - else
289   - % Select randomly new initial guesses within the specified
290   - % uncertainties.
291   - hfov0 = (hfovGuess - dhfov) + 2*dhfov*rand(1);
292   - lambda0 = (lambdaGuess - dlambda) + 2*dlambda*rand(1);
293   - phi0 = (phiGuess - dphi) + 2*dphi*rand(1);
294   - H0 = (HGuess - dH) + 2*dH*rand(1);
295   - theta0 = (thetaGuess - dtheta) + 2*dtheta*rand(1);
296   - end
297   -
298   - % Create vector cv0 for the initial guesses.
299   - i = 0;
300   - if dhfov > 0.0
301   - i = i+1;
302   - cv0(i) = hfov0;
303   - theOrder(i) = 1;
304   - end
305   - if dlambda > 0.0
306   - i = i + 1;
307   - cv0(i) = lambda0;
308   - theOrder(i) = 2;
309   - end
310   - if dphi > 0.0
311   - i = i + 1;
312   - cv0(i) = phi0;
313   - theOrder(i) = 3;
314   - end
315   - if dH > 0.0
316   - i = i + 1;
317   - cv0(i) = H0;
318   - theOrder(i) = 4;
319   - end
320   - if dtheta > 0.0
321   - i = i + 1;
322   - cv0(i) = theta0;
323   - theOrder(i) = 5;
324   - end
325   -
326   - [cv, errGeoFit] = fminsearch('g_error_geofit',cv0,options, ...
327   - imgWidth,imgHeight,xp,yp,ic,jc,...
328   - hfov,lambda,phi,H,theta,...
329   - hfov0,lambda0,phi0,H0,theta0,...
330   - hfovGuess,lambdaGuess,phiGuess,HGuess,thetaGuess,...
331   - dhfov,dlambda,dphi,dH,dtheta,...
332   - LON0,LAT0,...
333   - i_gcp,j_gcp,lon_gcp0,lat_gcp0,h_gcp,...
334   - theOrder,field);
335   -
336   - if errGeoFit < bestErrGeoFit
337   - bestErrGeoFit = errGeoFit;
338   - cvBest = cv;
339   - end
340   -
341   - fprintf('\n')
342   - fprintf(' Iteration # (iMinimize): %i\n',iMinimize);
343   - fprintf(' Max. number of iteration (nMinimize): %i\n',nMinimize);
344   - fprintf(' RSM error (m) for this iteration (errGeoFit): %f\n',errGeoFit);
345   - fprintf(' Best RSM error (m) so far (bestErrGeoFit): %f\n',bestErrGeoFit);
346   -
347   - end
348   -
349   - for i = 1:length(theOrder)
350   - if theOrder(i) == 1; hfov = cvBest(i); end
351   - if theOrder(i) == 2; lambda = cvBest(i); end
352   - if theOrder(i) == 3; phi = cvBest(i); end
353   - if theOrder(i) == 4; H = cvBest(i); end
354   - if theOrder(i) == 5; theta = cvBest(i); end
355   - end
356   -
357   - fprintf('\n')
358   - fprintf('PARAMETERS AFTER GEOMETRICAL RECTIFICATION \n')
359   - fprintf(' Field of view (hfov): %f\n',hfov)
360   - fprintf(' Dip angle (lambda): %f\n',lambda)
361   - fprintf(' Tilt angle (phi): %f\n',phi)
362   - fprintf(' Camera altitude (H): %f\n',H)
363   - fprintf(' View angle from North (theta): %f\n',theta)
364   - fprintf('\n')
365   -
366   - if length(theOrder) > 1
367   - fprintf('The rms error after geometrical correction (m): %f\n',bestErrGeoFit);
368   - end
369   -end
370   -
371   -% Project the GCP that have elevation.
372   -[lon_gcp,lat_gcp] = g_proj_GCP(LON0,LAT0,H,lon_gcp0,lat_gcp0,h_gcp,field);
373   -
374   -%%
375   -
376   -% Now construct the matrices LON and LAT for the entire image using the
377   -% camera parameters found by minimization just above.
378   -
379   -% Camera coordinate of all pixels
380   -xp = repmat([1:imgWidth]',1,imgHeight);
381   -yp = repmat([1:imgHeight],imgWidth,1);
382   -
383   -% Transform camera coordinate to ground coordinate.
384   -[LON LAT] = g_pix2ll(xp,yp,imgWidth,imgHeight,ic,jc,...
385   - hfov,lambda,phi,theta,H,LON0,LAT0,field);
386   -
387   -
388   -%% Apply polynomial correction if requested.
389   -if polyCorrection == true
390   - [LON LAT errPolyFit] = g_poly(LON,LAT,LON0,LAT0,i_gcp,j_gcp,lon_gcp,lat_gcp,polyOrder,field);
391   - fprintf('The rms error after polynomial stretching (m): %f\n',errPolyFit)
392   -else
393   - errPolyFit = NaN;
394   -end
395   -%%
396   -
397   -% Compute the effective resolution
398   -Delta = g_res(LON, LAT, field);
399   -
400   -fprintf('\n')
401   -fprintf('Saving rectification file in: %s\n',outputFname);
402   -
403   -save(outputFname,'imgFname','firstImgFname','lastImgFname',...
404   - 'LON','LAT',...
405   - 'LON0','LAT0',...
406   - 'lon_gcp0','lat_gcp0',...
407   - 'lon_gcp','lat_gcp','h_gcp',...
408   - 'i_gcp','j_gcp',...
409   - 'hfov','lambda','phi','H','theta',...
410   - 'errGeoFit','errPolyFit',...
411   - 'precision','Delta');
412   -
413   -clear LON LAT
414   -
415   -fprintf('\n')
416   -fprintf('Making figure\n');
417   -
418   -if field
419   - g_viz_field(imgFname,outputFname);
420   -else
421   - g_viz_lab(imgFname,outputFname);
422   -end
423   -
424   -print -dpng image2.png
  1 +function g_rect
  2 +%G_RECT - Main function for georectifying oblique images.
  3 +%
  4 +% Syntax: Simply type g_rect at the command line and follow the
  5 +% instructions.
  6 +%
  7 +% Inputs:
  8 +% The function G_RECT reads an input parameter file that contains all
  9 +% information required to perfrorm a georectification of a given image.
  10 +% The format of this parameter file is detailed on-line on the G_RECT
  11 +% Wiki page.
  12 +%
  13 +% Outputs:
  14 +% The function G_RECT creates an output file that contains the following
  15 +% variables:
  16 +%
  17 +% imgFname: The reference image for the georectification.
  18 +%
  19 +% firstImgFname: The first image of a sequence of images to which the
  20 +% georectification could be applied to. This is really
  21 +% just a comment.
  22 +%
  23 +% lastImgFname: The last image of a sequence of images to which the
  24 +% georectification could be applied to. This is really
  25 +% just a comment.
  26 +%
  27 +% LON: The main matrix that contain the longitude of each
  28 +% pixel of the referecne image (imgFname). Note that
  29 +% if the package is used to rectify lab images this
  30 +% matrix rather contains the distance in meters from
  31 +% a predefined x-origin.
  32 +%
  33 +% LAT: Same as LON but for the latitude.
  34 +%
  35 +% LON0: A scalar for the longitude of the camera or, in the
  36 +% case of a lab setup its cartesian coordinate in meters.
  37 +%
  38 +% LAT0: Same as LON0 but for the latitude.
  39 +%
  40 +% lon0_gcp: A vector containing the longitude of each ground
  41 +% control points (GCP). For the lab case this is the
  42 +% cartesian coordinate in meters. These GCPs may have
  43 +% elevations.
  44 +%
  45 +% lat0_gcp: Same as lon_gcp for latitude.
  46 +%
  47 +% lon_gcp: A vector containing the longitude of each ground
  48 +% control points (GCP) projected onto the water level
  49 +% i.e. at 0 m of elevation. For the lab case this is
  50 +% the cartesian coordinate in meters.
  51 +%
  52 +% lat_gcp: Same as lon_gcp for latitude.
  53 +%
  54 +% h_gcp: The elevation in meter of the GCPs. The elevation is
  55 +% 0 m if taken at water level.
  56 +%
  57 +% i_gcp: The horizontal index of the image ground control
  58 +% points.
  59 +%
  60 +% j_gcp: The vertical index of the image ground control
  61 +% points.
  62 +%
  63 +% hfov: The camera horizontal field of view [degree].
  64 +%
  65 +% phi: Camera tilt angle [degree].
  66 +%
  67 +% lambda: Camera dip angle [degree].
  68 +%
  69 +% H: The camera altitude relative to the water [m].
  70 +%
  71 +% theta: View angle clockwise from North [degree].
  72 +%
  73 +% errGeoFit: The rms error of the georectified image after
  74 +% geometrical transformation [m].
  75 +%
  76 +% errPolyFit: The rms error of the georectified image after
  77 +% geometrical transformation and the polynomial
  78 +% correction [m].
  79 +%
  80 +% precision: Calculation can be done in single or double
  81 +% precisions as defined by the user in the parameter
  82 +% file. With today's computers this is now obsolete
  83 +% and calculations can always be done in double
  84 +% precision.
  85 +%
  86 +% Delta: The effective resolution of the georectified image.
  87 +%
  88 +% Other m-files required: Works best with the m_map package for
  89 +% visualization.
  90 +% Subfunctions: all functions contained within the G_RECT folder.
  91 +%
  92 +% Author: Daniel Bourgault
  93 +% Institut des sciences de la mer de Rimouski
  94 +% email: daniel_bourgault@uqar.ca
  95 +% Website: https://srwebpolr01.uqar.ca/g_rect/
  96 +% February 2013
  97 +%
  98 +
  99 +%%
  100 +% The minimization is repeated nMinimize times where each time a random
  101 +% combination of the initial guesses is chosen within the given
  102 +% uncertainties provided by the user. This is becasue the algorithm often
  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.
  105 +nMinimize = 50;
  106 +
  107 +%% Read the parameter file
  108 +% Count the number of header lines before the ground control points (GCPs)
  109 +% The GCPs start right after the variable gcpData is set to true.
  110 +display(' ');
  111 +display(' Welcome to g_rect: a package for georectifying oblique images on a flat ocean');
  112 +display(' Authors: Daniel Bourgault and Rich Pawlowicz');
  113 +display(' ');
  114 +inputFname = input(' Enter the filename for the input parameters: ','s');
  115 +fid = fopen(inputFname);
  116 +
  117 +nHeaderLine = 0;
  118 +gcpData = false;
  119 +
  120 +% Read and execute each line of the parameter file until gcpData = true
  121 +% after which the GCP data appear and are read below with the importdata
  122 +% function.
  123 +while gcpData == false
  124 + eval(fgetl(fid));
  125 + nHeaderLine = nHeaderLine + 1;
  126 +end
  127 +fclose(fid);
  128 +
  129 +%% Import the GCP data at the end of the parameter file
  130 +gcp = importdata(inputFname,' ',nHeaderLine);
  131 +i_gcp = gcp.data(:,1);
  132 +j_gcp = gcp.data(:,2);
  133 +lon_gcp0 = gcp.data(:,3);
  134 +lat_gcp0 = gcp.data(:,4);
  135 +
  136 +[ngcp n_column] = size(gcp.data)
  137 +
  138 +% If there are 5 columns it means that elevation are provided
  139 +if n_column == 5
  140 + h_gcp = gcp.data(:,5);
  141 +else % otherwise elevation are considered to be zero
  142 + h_gcp(1:ngcp) = 0.0;
  143 +end
  144 +
  145 +%%
  146 +% Check if the elevation of the GCPs are not too high and above
  147 +% a certain fraction (gamma) of the camera height. If so, stop.
  148 +gamma = 0.75;
  149 +i_bad = find(h_gcp > gamma*(H+dH));
  150 +if length(i_bad) > 0
  151 + display([' ']);
  152 + display([' WARNING:']);
  153 + for i = 1:length(i_bad)
  154 + display([' The elevation of GCP #',num2str(i_bad(i)),' is greater than ',num2str(gamma),'*(H+dH).']);
  155 + end
  156 + display([' FIX AND RERUN.']);
  157 + return
  158 +end
  159 +
  160 +% Get the image size
  161 +imgInfo = imfinfo(imgFname);
  162 +imgWidth = imgInfo.Width;
  163 +imgHeight = imgInfo.Height;
  164 +
  165 +if precision == 'single'
  166 + imgWidth = single(imgWidth);
  167 + imgHeight = single(imgHeight);
  168 +end
  169 +
  170 +%% Display information
  171 +fprintf('\n')
  172 +fprintf(' INPUT PARAMETERS\n')
  173 +fprintf(' Image filename: (imgFname):........... %s\n',imgFname)
  174 +fprintf(' First image: (firstImgFname):......... %s\n',firstImgFname)
  175 +fprintf(' Last image: (lastImgFname):........... %s\n',lastImgFname)
  176 +fprintf(' Output filename: (outputFname):....... %s\n',outputFname);
  177 +fprintf(' Image width (imgWidth):............... %i\n',imgWidth)
  178 +fprintf(' Image width (imgHeight):.............. %i\n',imgHeight)
  179 +fprintf(' Camera longitude (LON0):.............. %f\n',LON0)
  180 +fprintf(' Camera latitude (LAT0):............... %f\n',LAT0)
  181 +fprintf(' Principal point offset (ic):.......... %f\n',ic)
  182 +fprintf(' Principal point offset (jc):.......... %f\n',jc)
  183 +fprintf(' Field of view (hfov):................. %f\n',hfov)
  184 +fprintf(' Dip angle (lambda):................... %f\n',lambda)
  185 +fprintf(' Tilt angle (phi):..................... %f\n',phi)
  186 +fprintf(' Camera altitude (H):.................. %f\n',H)
  187 +fprintf(' View angle from North (theta):........ %f\n',theta)
  188 +fprintf(' Uncertainty in hfov (dhfov):.......... %f\n',dhfov)
  189 +fprintf(' Uncertainty in dip angle (dlambda):... %f\n',dlambda)
  190 +fprintf(' Uncertainty in tilt angle (dphi):..... %f\n',dphi)
  191 +fprintf(' Uncertainty in altitude (dH):......... %f\n',dH)
  192 +fprintf(' Uncertainty in view angle (dtheta):... %f\n',dtheta)
  193 +fprintf(' Polynomial order (polyOrder):......... %i\n',polyOrder)
  194 +fprintf(' Number of GCPs (ngcp):................ %i\n',ngcp)
  195 +fprintf(' Precision (precision):................ %s\n',precision)
  196 +fprintf(' Field or lab (field=true; lab=false):. %i\n',field)
  197 +fprintf('\n')
  198 +
  199 +% Display the image with GCPs;
  200 +%image(imread(imgFname));
  201 +imagesc(imread(imgFname));
  202 +colormap(gray);
  203 +hold on
  204 +for i = 1:ngcp
  205 + plot(i_gcp(i),j_gcp(i),'r.');
  206 + text(i_gcp(i),j_gcp(i),[' ',num2str(i),'(',num2str(h_gcp(i)),')'],...
  207 + 'color','r',...
  208 + 'horizontalalignment','left',...
  209 + 'fontsize',10);
  210 +end
  211 +title('Ground Control Points Number (elevation in meters)','color','r');
  212 +daspect([1 1 1]);
  213 +print('-dpng','-r300',[imgFname(1:end-4),'_GCP.png']);
  214 +
  215 +fprintf('\n')
  216 +ok = input('Is it ok to proceed with the rectification (y/n): ','s');
  217 +if ok ~= 'y'
  218 + return
  219 +end
  220 +
  221 +%%
  222 +nUnknown = 0;
  223 +if dhfov > 0.0; nUnknown = nUnknown+1; end
  224 +if dlambda > 0.0; nUnknown = nUnknown+1; end
  225 +if dphi > 0.0; nUnknown = nUnknown+1; end
  226 +if dH > 0.0; nUnknown = nUnknown+1; end
  227 +if dtheta > 0.0; nUnknown = nUnknown+1; end
  228 +
  229 +if nUnknown > ngcp
  230 + fprintf('\n')
  231 + fprintf('WARNING: \n');
  232 + fprintf(' The number of GCPs is < number of unknown parameters.\n');
  233 + fprintf(' Program stopped.\n');
  234 + return
  235 +end
  236 +
  237 +% Check for consistencies between number of GCPs and order of the polynomial
  238 +% correction
  239 +ngcp = length(i_gcp);
  240 +if ngcp < 3*polyOrder
  241 + fprintf('\n')
  242 + fprintf('WARNING: \n');
  243 + fprintf(' The number of GCPs is inconsistent with the order of the polynomial correction.\n');
  244 + fprintf(' ngcp should be >= 3*polyOrder.\n');
  245 + fprintf(' Polynomial correction will not be applied.\n');
  246 + polyCorrection = false;
  247 +else
  248 + polyCorrection = true;
  249 +end
  250 +if polyOrder == 0
  251 + polyCorrection = false;
  252 +end
  253 +
  254 +%% This is the main section for the minimization algorithm
  255 +
  256 +if nUnknown > 0
  257 +
  258 + % Options for the fminsearch function. May be needed for some particular
  259 + % problems but in general the default values should work fine.
  260 + %options=optimset('MaxFunEvals',100000,'MaxIter',100000);
  261 + %options=optimset('MaxFunEvals',100000,'MaxIter',100000,'TolX',1.d-12,'TolFun',1.d-12);
  262 + options = [];
  263 +
  264 + % Only feed the minimization algorithm with the GCPs. xp and yp are the
  265 + % image coordinate of these GCPs.
  266 + xp = i_gcp;
  267 + yp = j_gcp;
  268 +
  269 + % This is the call to the minimization
  270 + bestErrGeoFit = Inf;
  271 +
  272 + % Save inital guesses in new variables.
  273 + hfovGuess = hfov;
  274 + lambdaGuess = lambda;
  275 + phiGuess = phi;
  276 + HGuess = H;
  277 + thetaGuess = theta;
  278 +
  279 + for iMinimize = 1:nMinimize
  280 +
  281 + % First guesses for the minimization
  282 + if iMinimize == 1
  283 + hfov0 = hfov;
  284 + lambda0 = lambda;
  285 + phi0 = phi;
  286 + H0 = H;
  287 + theta0 = theta;
  288 + else
  289 + % Select randomly new initial guesses within the specified
  290 + % uncertainties.
  291 + hfov0 = (hfovGuess - dhfov) + 2*dhfov*rand(1);
  292 + lambda0 = (lambdaGuess - dlambda) + 2*dlambda*rand(1);
  293 + phi0 = (phiGuess - dphi) + 2*dphi*rand(1);
  294 + H0 = (HGuess - dH) + 2*dH*rand(1);
  295 + theta0 = (thetaGuess - dtheta) + 2*dtheta*rand(1);
  296 + end
  297 +
  298 + % Create vector cv0 for the initial guesses.
  299 + i = 0;
  300 + if dhfov > 0.0
  301 + i = i+1;
  302 + cv0(i) = hfov0;
  303 + theOrder(i) = 1;
  304 + end
  305 + if dlambda > 0.0
  306 + i = i + 1;
  307 + cv0(i) = lambda0;
  308 + theOrder(i) = 2;
  309 + end
  310 + if dphi > 0.0
  311 + i = i + 1;
  312 + cv0(i) = phi0;
  313 + theOrder(i) = 3;
  314 + end
  315 + if dH > 0.0
  316 + i = i + 1;
  317 + cv0(i) = H0;
  318 + theOrder(i) = 4;
  319 + end
  320 + if dtheta > 0.0
  321 + i = i + 1;
  322 + cv0(i) = theta0;
  323 + theOrder(i) = 5;
  324 + end
  325 +
  326 + [cv, errGeoFit] = fminsearch('g_error_geofit',cv0,options, ...
  327 + imgWidth,imgHeight,xp,yp,ic,jc,...
  328 + hfov,lambda,phi,H,theta,...
  329 + hfov0,lambda0,phi0,H0,theta0,...
  330 + hfovGuess,lambdaGuess,phiGuess,HGuess,thetaGuess,...
  331 + dhfov,dlambda,dphi,dH,dtheta,...
  332 + LON0,LAT0,...
  333 + i_gcp,j_gcp,lon_gcp0,lat_gcp0,h_gcp,...
  334 + theOrder,field);
  335 +
  336 + if errGeoFit < bestErrGeoFit
  337 + bestErrGeoFit = errGeoFit;
  338 + cvBest = cv;
  339 + end
  340 +
  341 + fprintf('\n')
  342 + fprintf(' Iteration # (iMinimize): %i\n',iMinimize);
  343 + fprintf(' Max. number of iteration (nMinimize): %i\n',nMinimize);
  344 + fprintf(' RSM error (m) for this iteration (errGeoFit): %f\n',errGeoFit);
  345 + fprintf(' Best RSM error (m) so far (bestErrGeoFit): %f\n',bestErrGeoFit);
  346 +
  347 + end
  348 +
  349 + for i = 1:length(theOrder)
  350 + if theOrder(i) == 1; hfov = cvBest(i); end
  351 + if theOrder(i) == 2; lambda = cvBest(i); end
  352 + if theOrder(i) == 3; phi = cvBest(i); end
  353 + if theOrder(i) == 4; H = cvBest(i); end
  354 + if theOrder(i) == 5; theta = cvBest(i); end
  355 + end
  356 +
  357 + fprintf('\n')
  358 + fprintf('PARAMETERS AFTER GEOMETRICAL RECTIFICATION \n')
  359 + fprintf(' Field of view (hfov): %f\n',hfov)
  360 + fprintf(' Dip angle (lambda): %f\n',lambda)
  361 + fprintf(' Tilt angle (phi): %f\n',phi)
  362 + fprintf(' Camera altitude (H): %f\n',H)
  363 + fprintf(' View angle from North (theta): %f\n',theta)
  364 + fprintf('\n')
  365 +
  366 + if length(theOrder) > 1
  367 + fprintf('The rms error after geometrical correction (m): %f\n',bestErrGeoFit);
  368 + end
  369 +end
  370 +
  371 +% Project the GCP that have elevation.
  372 +[lon_gcp,lat_gcp] = g_proj_GCP(LON0,LAT0,H,lon_gcp0,lat_gcp0,h_gcp,field);
  373 +
  374 +%%
  375 +
  376 +% Now construct the matrices LON and LAT for the entire image using the
  377 +% camera parameters found by minimization just above.
  378 +
  379 +% Camera coordinate of all pixels
  380 +xp = repmat([1:imgWidth]',1,imgHeight);
  381 +yp = repmat([1:imgHeight],imgWidth,1);
  382 +
  383 +% Transform camera coordinate to ground coordinate.
  384 +[LON LAT] = g_pix2ll(xp,yp,imgWidth,imgHeight,ic,jc,...
  385 + hfov,lambda,phi,theta,H,LON0,LAT0,field);
  386 +
  387 +
  388 +%% Apply polynomial correction if requested.
  389 +if polyCorrection == true
  390 + [LON LAT errPolyFit] = g_poly(LON,LAT,LON0,LAT0,i_gcp,j_gcp,lon_gcp,lat_gcp,polyOrder,field);
  391 + fprintf('The rms error after polynomial stretching (m): %f\n',errPolyFit)
  392 +else
  393 + errPolyFit = NaN;
  394 +end
  395 +%%
  396 +
  397 +% Compute the effective resolution
  398 +Delta = g_res(LON, LAT, field);
  399 +
  400 +fprintf('\n')
  401 +fprintf('Saving rectification file in: %s\n',outputFname);
  402 +
  403 +save(outputFname,'imgFname','firstImgFname','lastImgFname',...
  404 + 'LON','LAT',...
  405 + 'LON0','LAT0',...
  406 + 'lon_gcp0','lat_gcp0',...
  407 + 'lon_gcp','lat_gcp','h_gcp',...
  408 + 'i_gcp','j_gcp',...
  409 + 'hfov','lambda','phi','H','theta',...
  410 + 'errGeoFit','errPolyFit',...
  411 + 'precision','Delta');
  412 +
  413 +clear LON LAT
  414 +
  415 +fprintf('\n')
  416 +fprintf('Making figure\n');
  417 +
  418 +if field
  419 + g_viz_field(imgFname,outputFname);
  420 +else
  421 + g_viz_lab(imgFname,outputFname);
  422 +end
  423 +print('-dpng','-r300',[imgFname(1:end-4),'_grect.png']);
... ...