Commit 56902974 authored by Daniel Bourgault's avatar Daniel Bourgault

Important changes made to the now obsolete variable field that has been...

Important changes made to the now obsolete variable field that has been changed to refFrame which is clearer. refFrame can now take the vamue Geodetic or Cartesian. Also much more comments added to the parameters.dat files. Finally, addition of Example_2 provided by Tom Shand.
parent c1a5bf26
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The image found in this folder ('IMG_6614.JPG') was taken from the ski resort
% Massif de Charlevoix in Quebec (Canada). It shows sea-ice in the St. Lawrence
% Estuary. This image is not associated with any publication and is provided here
% only as an example for the use of the g_rect package.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INPUT/OUTPUT INFORMATION
% imgFname: This is the only image that will actually be treated
%
% firstImgFname: This is really just a comment not actually used by the algorithm
% that indicates the name of the first image of a sequence to which
% the georectification of image 'imgFname' could also be applied to.
%
% lastImgFname: Same as 'firstImgFname'but for the last sequential image.
%
% outputFname: This is the name of the output fil that will contain, among other
% variables, the matrices LON and LAT that will give the ground
% coordinate of every pixel of the image that are below the horizon.
%
imgFname = 'IMG_6614.JPG';
firstImgFname = 'IMG_6614.JPG';
lastImgFname = 'IMG_6614.JPG';
outputFname = 'IMG_6614_grect.mat';
%% FRAME OF REFERENCE
% The frame of reference could be either 'Geodetic' or 'Cartesian'.
% If 'Geodetic' is used, longitudes and latitudes are expected for positions.
% If 'Cartesian' is used, x-y positions are expected (in m)
%
frameRef = 'Geodetic';
%% CAMERA POSITION
% Set the longitude and latitude for ‘Geodetic’ frame of reference
% Set the x-y coordinate for ‘Cartesian’ frame of reference
%
LON0 = -70.6010167;
LAT0 = 47.2713000;
%% OFFSET
% Offset from center of the principal point (generally zero)
%
ic = 0;
jc = 0;
%% CAMERA PARAMETERS
% These are either the initial guesses for the uncertain parameters, i.e. those
% with an uncertainty > 0 (see next section) or the exact values for the parameters
% with an uncertainty set to 0 below.
%
% hfov: Field of view of the camera (degree)
%
% lambda: Dip angle (degree) below horizontal (e.g. straight down = 90, horizontal = 0)
%
% phi: Tilt angle (degree), generally close to 0.
%
% theta: View angle (degree) clockwise from North (e.g. straight East = 90)
%
% H: Camera altitude relative to the local water level (m)
%
hfov = 65;
lambda = 2;
phi = 0;
theta = 70;
H = 720;
%% UNCERTAINTIES
% Set here the uncertainties of the associated camera parameters.
% Set the uncertainty to 0 for fixed parameters.
%
dhfov = 20;
dlambda = 10;
dphi = 5;
dtheta = 20;
dH = 0;
%% POLYNOMIAL CORRECTION
% After the geometrical correction, a polynomial correction of degree 1 or 2
% could be applied. This could correct for some unknown distortions that cannot be
% corrected on geometrical grounds. Play carefully with this option as it is a purely
% mathematical fit which may lead to unphysical corrections or may hide other
% prior problems with the actual geometrical correction. Always first set this option
% to '0' in which case there will be no polynomial correction applied. In principle,
% the geometrical fit without this option should already be pretty good. You could then
% fine tune the image with this option. Be extra careful when using the second order
% polynomial fit, especially outside the region of the GCPs as the image could there
% be completely distorted.
%
polyOrder = 0;
%% PRECISION
% To save memory, calculations can be done in single precision.
% For higher precision set the variable 'precision' to 'double';
%
precision = 'double';
%% GROUND CONTROL POINTS (GCPs).
% The GCP data must come right after the 'gcpData = true' instruction
% Column 1: horizontal image index of GCPs
% Column 2: vertical image index of GCPs
% Column 3: longitude (Geodetic) or x-position (Cartesian) of GCPs
% Column 4: latitude (Geodetic) or y-position (Cartesian) of GCPs
% Column 5 (optional): elevation (in m) of GCPs.
gcpData = true;
360 829 -70.561367 47.303783 0
54 719 -70.54500 47.335 0
99 661 -70.505 47.375 0
452 641 -70.435 47.389 0
429 633 -70.418 47.408 0
816 644 -70.393 47.368 0
\ No newline at end of file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% The image found in this folder ('Camera01.jpg') as well as the data found below in
% this parameter file were graciously provided by Tom Shand from Tonkin + Taylor Ltd.
% Auckland, New Zealand (www.tonkintaylor.co.nz).
%
% The reference for the image and the data is:
%
% Shand, T., Reinen-Hamill, R., Weppe, S. and A. Short (2019) Development of a
% framework for assessing effects of coastal engineering works on a surf break.
% Australasian Coasts and Ports Conference, September 2019, Hobart, Australia.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%% INPUT/OUTPUT INFORMATION
% imgFname: This is the only image that will actually be treated
%
% firstImgFname: This is really just a comment not actually used by the algorithm
% that indicates the name of the first image of a sequence to which
% the georectification of image 'imgFname' could also be applied to.
%
% lastImgFname: Same as 'firstImgFname'but for the last sequential image.
%
% outputFname: This is the name of the output fil that will contain, among other
% variables, the matrices LON and LAT that will give the ground
% coordinate of every pixel of the image that are below the horizon.
%
imgFname = 'Camera01.jpg';
firstImgFname = 'Camera01.jpg';
lastImgFname = 'Camera01.jpg';
outputFname = 'Camera01_grect.mat';
%% FRAME OF REFERENCE
% The frame of reference could be either 'Geodetic' or 'Cartesian'.
% If 'Geodetic' is used, longitudes and latitudes are expected for positions.
% If 'Cartesian' is used, x-y positions are expected (in m)
%
frameRef = 'Cartesian';
%% CAMERA POSITION
% Set the longitude and latitude for ‘Geodetic’ frame of reference
% Set the x-y coordinate for ‘Cartesian’ frame of reference.
%
LON0 = 1661768;
LAT0 = 5316087;
%% OFFSET
% Offset from center of the principal point (generally zero)
%
ic = 0;
jc = 0;
%% CAMERA PARAMETERS
% These are either the initial guesses for the uncertain parameters, i.e. those
% with an uncertainty > 0 (see next section) or the exact values for the parameters
% with an uncertainty set to 0 below.
%
% hfov: Field of view of the camera (degree)
%
% lambda: Dip angle (degree) below horizontal (e.g. straight down = 90, horizontal = 0)
%
% phi: Tilt angle (degree), generally close to 0.
%
% theta: View angle (degree) clockwise from North (e.g. straight East = 90)
%
% H: Camera altitude relative to the local water level (m)
%
hfov = 70;
lambda = 20;
phi = 0;
theta = 90;
H = 30;
%% UNCERTAINTIES
% Set here the uncertainties of the associated camera parameters.
% Set the uncertainty to 0 for fixed parameters.
%
dhfov = 10;
dlambda = 20;
dphi = 2;
dtheta = 10;
dH = 0;
%% POLYNOMIAL CORRECTION
% After the geometrical correction, a polynomial correction of degree 1 or 2
% could be applied. This could correct for some unknown distortions that cannot be
% corrected on geometrical grounds. Play carefully with this option as it is a purely
% mathematical fit which may lead to unphysical corrections or may hide other
% prior problems with the actual geometrical correction. Always first set this option
% to '0' in which case there will be no polynomial correction applied. In principle,
% the geometrical fit without this option should already be pretty good. You could then
% fine tune the image with this option. Be extra careful when using the second order
% polynomial fit, especially outside the region of the GCPs as the image could there
% be completely distorted.
%
polyOrder = 0;
%% PRECISION
% To save memory, calculations can be done in single precision.
% For higher precision set the variable 'precision' to 'double';
%
precision = 'double';
%% GROUND CONTROL POINTS (GCPs).
% The GCP data must come right after the 'gcpData = true' instruction
% Column 1: horizontal image index of GCPs
% Column 2: vertical image index of GCPs
% Column 3: longitude (Geodetic) or x-position (Cartesian) of GCPs
% Column 4: latitude (Geodetic) or y-position (Cartesian) of GCPs
% Column 5 (optional): elevation (in m) of GCPs.
%
gcpData = true;
1818 490 1661903.38 5316012.52 6.28
1359 662 1661846.11 5316069.20 6.11
543 1080 1661802.48 5316102.36 6.27
1584 479 1661926.61 5316026.49 3.90
1226 535 1661902.92 5316070.64 2.45
617 836 1661828.83 5316109.15 2.47
1475 347 1662214.16 5315956.24 1.58
974 355 1662195.27 5316115.86 1.10
14 583 1661893.51 5316195.60 0.43
595 382 1662132.48 5316216.56 0.60
132 408 1662085.97 5316320.83 0.42
\ No newline at end of file
......@@ -8,7 +8,7 @@ outputFname = 'IMG_8368_grect.mat';
% Field or lab case situation.
% Set field = true for field situation and field = false for lab situation.
field = false;
frameRef = 'Cartesian';
% Camera position
% lat/lon for field situation
......
% I/O information
imgFname = 'IMG_6614.JPG';
firstImgFname = 'IMG_6614.JPG';
lastImgFname = 'IMG_6614.JPG';
outputFname = 'IMG_6614_grect.mat';
% Field or lab case situation.
% Set field = true for field situation and field = false for lab situation.
field = true;
% Camera position
% lat/lon for field situation
% meter for lab situation
LON0 = -70.6010167;
LAT0 = 47.2713000;
% Offset from center of the principal point (generally zero)
ic = 0;
jc = 0;
% Parameters
hfov = 65.0; % Field of view of the camera
lambda = 2; % Dip angle above vertical (e.g. straight down = 90, horizontal = 0)
phi = 0.0; % Tilt angle (generally close to 0).
H = 720; % Camera altitude
theta = 70.0; % View angle clockwise from North (e.g. straight East = 90)
% Uncertainty in parameters. Set the uncertainty to 0.0 for fixed parameters.
dhfov = 20.0;
dlambda = 10.0;
dphi = 5.0;
dH = 0.0;
dtheta = 20.0;
% Order of the polynomial correction (0, 1 or 2)
polyOrder = 0;
% To save memory calculation can be done in single precision.
% For higher precision set the variable 'precision' to 'double';
precision = 'double';
% Ground Control Points (GCP).
% The data must come right after the gcpData = true
gcpData = true;
360 829 -70.561367 47.303783 0
54 719 -70.54500 47.335 0
99 661 -70.505 47.375 0
452 641 -70.435 47.389 0
429 633 -70.418 47.408 0
816 644 -70.393 47.368 0
\ No newline at end of file
function distance = g_dist(lon1,lat1,lon2,lat2,field);
function distance = g_dist(lon1,lat1,lon2,lat2,frameRef);
% This function computes the distance (m) between two points on a
% Cartesian Earth given their lat-lon coordinate.
%
% If field = false then simple cartesian transformation
%
% This function computes the distance (m) between two points given
% either their geodetic coordinates or cartesian coordinates.
%
dlon = lon2 - lon1;
dlat = lat2 - lat1;
if field
meterPerDegLat = 1852*60.0;
meterPerDegLon = meterPerDegLat * cosd(lat1);
dx = dlon*meterPerDegLon;
dy = dlat*meterPerDegLat;
else
if strcmp(frameRef,'Geodetic') == true
dx = dlon;
dy = dlat;
% Simplest calculation on a mercator plane
%dlon = lon2 - lon1;
%dlat = lat2 - lat1;
% meterPerDegLat = 1852*60.0;
% meterPerDegLon = meterPerDegLat * cosd(lat1);
% dx = dlon*meterPerDegLon;
% dy = dlat*meterPerDegLat;
% distance = sqrt(dx^2 + dy^2)
% Using the m_map package on a spherical Earth (more precise)
distance = m_lldist([lon2 lon1],[lat2, lat1])*1000;
end
distance = sqrt(dx^2 + dy^2);
elseif strcmp(frameRef,'Cartesian') == true
dx = lon2 - lon1;
dy = lat2 - lat1;
distance = sqrt(dx^2 + dy^2);
end
......@@ -5,10 +5,10 @@ function errGeoFit = g_error_geofit(cv,imgWidth,imgHeight,xp,yp,ic,jc,...
dhfov,dlambda,dphi,dH,dtheta,...
LON0,LAT0,...
i_gcp,j_gcp,lon_gcp0,lat_gcp0,h_gcp,...
theOrder,field);
theOrder,frameRef);
for i=1:length(theOrder)
for i = 1:length(theOrder)
if theOrder(i) == 1; hfov = cv(i); end
if theOrder(i) == 2; lambda = cv(i); end
if theOrder(i) == 3; phi = cv(i); end
......@@ -18,7 +18,7 @@ end
% Perform the geometrical transformation
[LON LAT] = g_pix2ll(xp,yp,imgWidth,imgHeight,ic,jc,...
hfov,lambda,phi,theta,H,LON0,LAT0,field);
hfov,lambda,phi,theta,H,LON0,LAT0,frameRef);
% Calculate the error between ground control points (GCPs)
% and image control points once georectified.
......@@ -35,12 +35,12 @@ ngcpFinite = 0;
errGeoFit = 0.0;
% Project the GCPs onto the water surface given their elevation
[lon_gcp,lat_gcp] = g_proj_GCP(LON0,LAT0,H,lon_gcp0,lat_gcp0,h_gcp,field);
[lon_gcp,lat_gcp] = g_proj_GCP(LON0,LAT0,H,lon_gcp0,lat_gcp0,h_gcp,frameRef);
for k = 1:ngcp
% Calculate the distance (i.e. error) between GCP and rectificed ICPs.
distance = g_dist(lon_gcp(k),lat_gcp(k),LON(k),LAT(k),field);
distance = g_dist(lon_gcp(k),lat_gcp(k),LON(k),LAT(k),frameRef);
% Check if the distance is finite. The distance may be NaN for some
% GCPs that may temporarily be above the horizon. Those points are
......@@ -64,4 +64,3 @@ if abs(lambda - lambdaGuess) > dlambda; errGeoFit = inf; end
if abs(phi - phiGuess) > dphi; errGeoFit = inf; end
if abs(H - HGuess) > dH; errGeoFit = inf; end
if abs(theta - thetaGuess) > dtheta; errGeoFit = inf; end
function err=g_error_polyfit(cv,lonlon,latlat,err_ll,order);
%function err=g_error_polyfit(cv,lonlon,latlat,err_ll,order);
function err = g_error_polyfit(cv,lonlon,latlat,err_ll,order);
n_gcp=length(err_ll);
err=0;
n_gcp = length(err_ll);
err = 0;
for k=1:n_gcp
for k = 1:n_gcp
if order == 1
......
function [LON,LAT] = g_pix2ll(xp,yp,imgWidth,imgHeight,ic,jc,...
hfov,lambda,phi,theta,H,LON0,LAT0,field);
hfov,lambda,phi,theta,H,LON0,LAT0,frameRef);
% G_PIX2LL Converts pixel to ground coordinates
%
% input:
......@@ -13,7 +13,7 @@ function [LON,LAT] = g_pix2ll(xp,yp,imgWidth,imgHeight,ic,jc,...
% phi: Tilt angle clockwise around the principal axis
% theta: View angle clockwise from North (e.g. East = 90)
% H: Camera altitude (m) above surface of interest.
% LON0, LAT0: Camera longitude and latitude position
% LON0, LAT0: Camera geodetic coordinates or cartesian coordinates
%
% output: LAT,LON: Ground coordinates
%
......@@ -123,10 +123,10 @@ x_w = x_w.*s2f;
z_w = z_w.*s2f;
% Convert coordinates to lat/lon using locally cartesian assumption
if field
if strcmp(frameRef,'Geodetic') == true
LON = x_w/meterPerDegLon + LON0;
LAT = z_w/meterPerDegLat + LAT0;
else
elseif strcmp(frameRef,'Cartesian') == true
LON = x_w + LON0;
LAT = z_w + LAT0;
end
function [LONpc, LATpc, err_polyfit] = g_poly(LON,LAT,LON0,LAT0,i_gcp,j_gcp,lon_gcp,lat_gcp,p_order,field);
function [LONpc, LATpc, err_polyfit] = g_poly(LON,LAT,LON0,LAT0,i_gcp,j_gcp,lon_gcp,lat_gcp,p_order,frameRef);
ngcp=length(i_gcp);
% LON0 and LAT0 are removed to facilitate the fiminsearch algorithm
for k = 1:ngcp
lonlon(k) = LON(i_gcp(k),j_gcp(k))-LON0;
latlat(k) = LAT(i_gcp(k),j_gcp(k))-LAT0;
err_lon(k) = lonlon(k)-(lon_gcp(k)-LON0);
err_lat(k) = latlat(k)-(lat_gcp(k)-LAT0);
Delta_lon = LON(i_gcp(k)+1,j_gcp(k)) - LON(i_gcp(k),j_gcp(k));
Delta_lat = LAT(i_gcp(k)+1,j_gcp(k)) - LAT(i_gcp(k),j_gcp(k));
lonlon(k) = LON(i_gcp(k),j_gcp(k))-LON0;
latlat(k) = LAT(i_gcp(k),j_gcp(k))-LAT0;
err_lon(k) = lonlon(k)-(lon_gcp(k)-LON0);
err_lat(k) = latlat(k)-(lat_gcp(k)-LAT0);
if (abs(err_lon) - Delta_lon) < 0
err_lon(k) = 0.0;
end
if (abs(err_lat) - Delta_lat) < 0
err_lat(k) = 0.0;
end
Delta_lon = LON(i_gcp(k)+1,j_gcp(k)) - LON(i_gcp(k),j_gcp(k));
Delta_lat = LAT(i_gcp(k)+1,j_gcp(k)) - LAT(i_gcp(k),j_gcp(k));
if (abs(err_lon) - Delta_lon) < 0
err_lon(k) = 0.0;
end
if (abs(err_lat) - Delta_lat) < 0
err_lat(k) = 0.0;
end
end
LON = LON-LON0; LAT = LAT-LAT0;
LON = LON - LON0;
LAT = LAT - LAT0;
options = optimset('MaxFunEvals',1000000,'MaxIter',1000000,'TolFun',1.d-12,'TolX',1.d-12);
......@@ -62,7 +64,8 @@ elseif p_order == 2
err_lon(k) = lonlon(k)-(lon_gcp(k)-LON0);
err_lat(k) = latlat(k)-(lat_gcp(k)-LAT0);
end
LON = LON-LON0; LAT = LAT-LAT0;
LON = LON - LON0;
LAT = LAT - LAT0;
a(1:6) = 0.0;
b(1:6) = 0.0;
......@@ -70,8 +73,8 @@ elseif p_order == 2
a = fminsearch('g_error_polyfit',a,options,lonlon,latlat,err_lon,2);
b = fminsearch('g_error_polyfit',b,options,lonlon,latlat,err_lat,2);
b = fminsearch('g_error_polyfit',b,options,lonlon,latlat,err_lat,2);
LONpc=LON-(a(1)*LON.^2+a(2)*LAT.^2+a(3)*LON.*LAT+a(4)*LON+a(5)*LAT+a(6));
LATpc=LAT-(b(1)*LON.^2+b(2)*LAT.^2+b(3)*LON.*LAT+b(4)*LON+b(5)*LAT+b(6));
LONpc = LON-(a(1)*LON.^2+a(2)*LAT.^2+a(3)*LON.*LAT+a(4)*LON+a(5)*LAT+a(6));
LATpc = LAT-(b(1)*LON.^2+b(2)*LAT.^2+b(3)*LON.*LAT+b(4)*LON+b(5)*LAT+b(6));
end
......@@ -83,9 +86,8 @@ LATpc = LATpc + LAT0;
err_polyfit=0;
for k = 1:ngcp
DX = g_dist(LONpc(i_gcp(k),j_gcp(k)),LATpc(i_gcp(k),j_gcp(k)),lon_gcp(k),lat_gcp(k),field);
DX = g_dist(LONpc(i_gcp(k),j_gcp(k)),LATpc(i_gcp(k),j_gcp(k)),lon_gcp(k),lat_gcp(k),frameRef);
err_polyfit = err_polyfit + DX^2;
end
err_polyfit = sqrt(err_polyfit/ngcp);
err_polyfit = sqrt(err_polyfit/ngcp);
\ No newline at end of file
function [lon_gcp,lat_gcp] = g_proj_GCP(LON0,LAT0,H,lon_gcp0,lat_gcp0,h_gcp,field)
function [lon_gcp,lat_gcp] = g_proj_GCP(LON0,LAT0,H,lon_gcp0,lat_gcp0,h_gcp,frameRef)
%G_PROJ_GCP
%
% This function projects the GCPs that have an elevation greater than 0
......@@ -13,23 +13,23 @@ function [lon_gcp,lat_gcp] = g_proj_GCP(LON0,LAT0,H,lon_gcp0,lat_gcp0,h_gcp,fiel
%
% INPUTS:
%
% LON0: The longitiude of the camera
% LAT0: The latitude of the camera
% LON0: The longitiude or x-position of the camera
% LAT0: The latitude or y-position of the camera
% H: The elevation of the camera
% lon_gcp0: The longitudes of the GCPs
% lat_gcp0: The latitudes of the GCPs
% lon_gcp0: The longitudes or x-position of the GCPs
% lat_gcp0: The latitudes or y-position of the GCPs
% h_gcp: The elevations of the GCPs
%
% OUTPUTS:
% lon_gcp: The projected longitudes of the GCPs
% lat_gcp: The projected latitudes of the GCPs
% lon_gcp: The projected longitudes or x-position of the GCPs
% lat_gcp: The projected latitudes or y-position of the GCPs
%
n_gcp = length(h_gcp);
for i = 1:n_gcp
for i = 1:n_gcp
if field == true
if strcmp(frameRef,'Geodetic') == true
% The calculation for the distance between the GCPs and the camera
% is done on an elliptical earth and this takes a long calculation
......@@ -50,9 +50,9 @@ for i = 1:n_gcp
lat_gcp(i) = lat_gcp0(i);
end
else
elseif strcmp(frameRef,'Cartesian') == true
% For lab situations where positions are provided in meters rather
% For cartesian situations where positions are provided in meters rather
% than in latitude longitude the calculation is simpler.
dx = lon_gcp0(i) - LON0;
dy = lat_gcp0(i) - LAT0;
......
......@@ -24,34 +24,32 @@ function g_rect
% 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.
% frameRef: The reference frame used. It could be either
% 'Geodetic' or 'Cartesian'.
%
% LAT: Same as LON but for the latitude.
% LON: The main matrix that contain either the longitude of each
% pixel of the reference image (imgFname) or the x-coordinate
% (in m) depending on whether the package is used in geodetic or
% cartesian coordinates.
%
% LON0: A scalar for the longitude of the camera or, in the
% case of a lab setup its cartesian coordinate in meters.
% LAT: Same as LON but for the latitude or y-position.
%
% LON0: A scalar for the longitude or x-position of the camera.
%
% 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.
% lon0_gcp: A vector containing the longitude or x-position of each
% ground control points (GCP).
%
% lat0_gcp: Same as lon_gcp for latitude.
% lat0_gcp: Same as lon_gcp for latitude or y-position.
%
% 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.
% lon_gcp: A vector containing the longitude or x-position of each
% ground control points (GCP) projected onto the water level
% i.e. at 0 m of elevation.
%
% lat_gcp: Same as lon_gcp for latitude.
% lat_gcp: Same as lon_gcp for latitude or y-position.
%
% h_gcp: The elevation in meter of the GCPs. The elevation is
% h_gcp: The elevation (in m) of the GCPs. The elevation is
% 0 m if taken at water level.
%
% i_gcp: The horizontal index of the image ground control
......@@ -83,18 +81,17 @@ function g_rect
% and calculations can always be done in double
% precision.
%
% Delta: The effective resolution of the georectified image.
% Delta: The effective resolution (im m) 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.
% Other m-files required: The m_map package.
% Subfunctions: all functions contained within the subdirectories g_rect/src/
%
% 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
%
% email: daniel_bourgault@uqar.ca
% Initial development: February 2013
% Important update: May 2020
%%
% The minimization is repeated nMinimize times where each time a random
......@@ -118,7 +115,7 @@ 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
% after which the GCP data appear and are read below with the 'importdata'
% function.
while gcpData == false
eval(fgetl(fid));
......@@ -126,6 +123,27 @@ while gcpData == false
end
fclose(fid);
% The older version of the code had a variable called 'field' that could be
% set to 'true' or 'false' in the input parameters file depending on whether
% the application was done for a field situation (field = true) or a lab
% situation (field = false). It was implicitly assumed that positions were
% given by a geodetic system (lon-lat) for field situations and by a
% cartesian system (x-y) for lab situations. But this was a little confusing
% as field situations could also use Cartesian coordinates, for example when
% the Universal Transverse Mercator (UTM) system is used. This new version
% of the code now rather specifies the frame of reference used that could
% be either 'Cartesian' or 'Geodetic'. The following lines of code are
% simply here to allow this new version of the code to be compatible with
% the older version of the parameters file in which the now obsolete
% variable 'field' appeared and the current variable 'frameRef' was absent.
if ~exist('frameRef')
if field == true
frameRef = 'Geodetic';
elseif field == false
frameRef = 'Cartesian';
end
end
%% Import the GCP data at the end of the parameter file
gcp = importdata(inputFname,' ',nHeaderLine);
i_gcp = gcp.data(:,1);
......@@ -177,8 +195,9 @@ 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(' Frame of reference:................... %s\n',frameRef)
fprintf(' Camera longitude or x coord. (LON0):.. %f\n',LON0)
fprintf(' Camera latitude or y coord. (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)
......@@ -194,7 +213,6 @@ 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')