diff --git a/Examples/Field/IMG_6614.JPG b/Examples/Example_1/IMG_6614.JPG similarity index 100% rename from Examples/Field/IMG_6614.JPG rename to Examples/Example_1/IMG_6614.JPG diff --git a/Examples/Example_1/parameters.dat b/Examples/Example_1/parameters.dat new file mode 100755 index 0000000000000000000000000000000000000000..00ac84b2c0b2d63a4bcd0b69b7a4bd946200f126 --- /dev/null +++ b/Examples/Example_1/parameters.dat @@ -0,0 +1,112 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% 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 diff --git a/Examples/Example_2/Camera01.jpg b/Examples/Example_2/Camera01.jpg new file mode 100644 index 0000000000000000000000000000000000000000..a1cc3e81dc1481422bdd750c032e2c4f301d600e Binary files /dev/null and b/Examples/Example_2/Camera01.jpg differ diff --git a/Examples/Example_2/parameters.dat b/Examples/Example_2/parameters.dat new file mode 100644 index 0000000000000000000000000000000000000000..671f54c6ada15d0b162afd4047494f40e2ce4680 --- /dev/null +++ b/Examples/Example_2/parameters.dat @@ -0,0 +1,123 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% 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 diff --git a/Examples/Lab/IMG_8368.jpg b/Examples/Example_3/IMG_8368.jpg similarity index 100% rename from Examples/Lab/IMG_8368.jpg rename to Examples/Example_3/IMG_8368.jpg diff --git a/Examples/Lab/parameters.dat b/Examples/Example_3/parameters.dat similarity index 94% rename from Examples/Lab/parameters.dat rename to Examples/Example_3/parameters.dat index adf3a58e385f022aa3dcb978f8d3de9cf34e749f..1af87733fe6a046cb7714ab6074c80a9a70b7494 100755 --- a/Examples/Lab/parameters.dat +++ b/Examples/Example_3/parameters.dat @@ -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 diff --git a/Examples/Field/parameters.dat b/Examples/Field/parameters.dat deleted file mode 100755 index fc15b2fa4a10bb9696bc24a942e13eb532a7c269..0000000000000000000000000000000000000000 --- a/Examples/Field/parameters.dat +++ /dev/null @@ -1,50 +0,0 @@ -% 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 diff --git a/src/g_rect/g_dist.m b/src/g_rect/g_dist.m index 8187ac05f7ecc65aa91092e89c806c56b7e1780c..f4680ccf6971f7f69a6338a7fdebbe877456a7a1 100755 --- a/src/g_rect/g_dist.m +++ b/src/g_rect/g_dist.m @@ -1,29 +1,28 @@ -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 diff --git a/src/g_rect/g_error_geofit.m b/src/g_rect/g_error_geofit.m index 5df99e218a4659d83d5b23ad229ae8743ef569a6..1da9f1f61b10f5bf5860a52d4d5c3d69d33a3933 100755 --- a/src/g_rect/g_error_geofit.m +++ b/src/g_rect/g_error_geofit.m @@ -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 - diff --git a/src/g_rect/g_error_polyfit.m b/src/g_rect/g_error_polyfit.m index 44ea37d97776edb76955310cfaf284c36f636929..28fc1780a48707faae79e2d0e6d9d43f365e3ccb 100755 --- a/src/g_rect/g_error_polyfit.m +++ b/src/g_rect/g_error_polyfit.m @@ -1,10 +1,9 @@ -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 diff --git a/src/g_rect/g_pix2ll.m b/src/g_rect/g_pix2ll.m index 8a26c3f7ce3dd17314c2fdc4a3a9815d31782a45..4102de89c185e92c82e4329b13efb5e2040de18f 100755 --- a/src/g_rect/g_pix2ll.m +++ b/src/g_rect/g_pix2ll.m @@ -1,5 +1,5 @@ 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 diff --git a/src/g_rect/g_poly.m b/src/g_rect/g_poly.m index 2c8e29d8d4b0da227210ad5374a72627a64c482c..6283483463a0cbdbe8f9557b7b5685365ed7ce31 100755 --- a/src/g_rect/g_poly.m +++ b/src/g_rect/g_poly.m @@ -1,26 +1,28 @@ -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 diff --git a/src/g_rect/g_proj_GCP.m b/src/g_rect/g_proj_GCP.m index 900809bf5fdd0d7c2512147b0869e6ff64de956b..0ef10c78b193c824d3e4d7789af085ee0898f18f 100644 --- a/src/g_rect/g_proj_GCP.m +++ b/src/g_rect/g_proj_GCP.m @@ -1,4 +1,4 @@ -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; diff --git a/src/g_rect/g_rect.m b/src/g_rect/g_rect.m index e376a108df1257b9efb572a17c03cb7a56417610..562bcd81fd1111bc3dc34c302f9a8d486ad5c89c 100644 --- a/src/g_rect/g_rect.m +++ b/src/g_rect/g_rect.m @@ -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') % Display the image with GCPs; @@ -334,7 +352,7 @@ if nUnknown > 0 dhfov,dlambda,dphi,dH,dtheta,... LON0,LAT0,... i_gcp,j_gcp,lon_gcp0,lat_gcp0,h_gcp,... - theOrder,field); + theOrder,frameRef); if errGeoFit < bestErrGeoFit bestErrGeoFit = errGeoFit; @@ -372,7 +390,7 @@ if nUnknown > 0 end % Project the GCP that have 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); %% @@ -385,12 +403,12 @@ 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); + hfov,lambda,phi,theta,H,LON0,LAT0,frameRef); %% 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); + [LON LAT errPolyFit] = g_poly(LON,LAT,LON0,LAT0,i_gcp,j_gcp,lon_gcp,lat_gcp,polyOrder,frameRef); fprintf('The rms error after polynomial stretching (m): %f\n',errPolyFit) else errPolyFit = NaN; @@ -398,12 +416,13 @@ end %% % Compute the effective resolution -Delta = g_res(LON, LAT, field); +Delta = g_res(LON, LAT, frameRef); fprintf('\n') fprintf('Saving rectification file in: %s\n',outputFname); save(outputFname,'imgFname','firstImgFname','lastImgFname',... + 'frameRef',... 'LON','LAT',... 'LON0','LAT0',... 'lon_gcp0','lat_gcp0',... @@ -418,9 +437,9 @@ clear LON LAT fprintf('\n') fprintf('Making figure\n'); -if field - g_viz_field(imgFname,outputFname); -else - g_viz_lab(imgFname,outputFname); +if strcmp(frameRef,'Geodetic') + g_viz_geodetic(imgFname,outputFname); +elseif strcmp(frameRef,'Cartesian') + g_viz_cartesian(imgFname,outputFname); end print('-dpng','-r300',[imgFname(1:end-4),'_grect.png']); diff --git a/src/g_rect/g_res.m b/src/g_rect/g_res.m index fcb01093c6faf8902e25c6112c7d387ccafd3a82..2a54781d2087eb78e6a3df85ce8a5700bf3000c0 100755 --- a/src/g_rect/g_res.m +++ b/src/g_rect/g_res.m @@ -1,32 +1,29 @@ -function res = g_res_field(LON, LAT, field) +function res = g_res_field(LON, LAT, frameRef) % g_res_field - Function calculating the loss of resolution with de distance % from the camera % % Inputs: % -% LON : Longitude matrix obtained from g_rect -% LAT : Latitude matrix obtained from g_rect -% field : 'true' for field case. 'false' for lab case. +% LON : Longitude or x-positon matrix obtained from g_rect +% LAT : Latitude or y-position matrix obtained from g_rect +% frameRef : 'Geodetic' of 'Cartesian'. % % Outputs: % % res : The field of resolution which is degrading with distance from % the camera % -% 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 -% -% Author: Elie Dumas-Lefebvre -% Institut des Science de la Mer de Rimouski -% -% Note: Matricial formulation of the calculation for a faster execution +% Authors: +% Daniel Bourgault +% Institut des sciences de la mer de Rimouski +% email: daniel_bourgault@uqar.ca +% February 2013 % -% email: elie.dumas-lefebvre@uqar.ca -% March 2019 +% Elie Dumas-Lefebvre +% Institut des Science de la Mer de Rimouski +% Note: Matricial formulation of the calculation for a faster execution +% email: elie.dumas-lefebvre@uqar.ca +% March 2019 % % LAT and LON difference in the x-axis @@ -55,7 +52,8 @@ meterPerDegLat = 1852*60.0; meterPerDegLoni = meterPerDegLat.*cosd(LAT_meani); meterPerDegLonj = meterPerDegLat.*cosd(LAT_meanj); -if field + +if strcmp(frameRef,'Geodetic') == true % Conversion from degrees to meters dxi = dLONi.*meterPerDegLoni; @@ -63,7 +61,7 @@ if field dyi = dLATi.*meterPerDegLat; dyj = dLATj.*meterPerDegLat; -else +elseif strcmp(frameRef,'Cartesian') == true % Lab case. Data are already in meters. No conversion required. dxi = dLONi; @@ -78,6 +76,4 @@ deltai = sqrt(dxi.^2 + dyi.^2); deltaj = sqrt(dxj.^2 + dyj.^2); % field of resolution -res = sqrt(deltai.^2 + deltaj.^2); - -end \ No newline at end of file +res = sqrt(deltai.^2 + deltaj.^2); \ No newline at end of file diff --git a/src/g_rect/g_viz_lab.m b/src/g_rect/g_viz_cartesian.m similarity index 89% rename from src/g_rect/g_viz_lab.m rename to src/g_rect/g_viz_cartesian.m index 55b45117520933bb42533973aed7bf431ef7fe66..a1f9744bab9decef72c08d2b23c2b22f9b6de63b 100644 --- a/src/g_rect/g_viz_lab.m +++ b/src/g_rect/g_viz_cartesian.m @@ -1,4 +1,4 @@ -function g_viz_lab(imgFname,rectFile); +function g_viz_cartesian(imgFname,rectFile); figure @@ -33,16 +33,18 @@ end clear rgb0; int = int'; -hold on; +%hold on; colormap(gray); pcolor(LON,LAT,int); shading('flat'); +hold; + plot(LON0,LAT0,'kx','markersize',ms,'linewidth',lw); % Camera location %% Plot GCPs and ICPs. -for n=1:length(i_gcp) +for n = 1:length(i_gcp) % Plot the original GCPs that may have elevations plot(lon_gcp0(n),lat_gcp0(n),'ko','markersize',ms,'linewidth',lw); @@ -52,6 +54,7 @@ for n=1:length(i_gcp) % Plot the Image Control Points once georectified plot(LON(i_gcp(n),j_gcp(n)),LAT(i_gcp(n),j_gcp(n)),'rx','MarkerSize',ms,'linewidth',lw); + end xlabel('x (m)'); diff --git a/src/g_rect/g_viz_field.m b/src/g_rect/g_viz_geodetic.m similarity index 94% rename from src/g_rect/g_viz_field.m rename to src/g_rect/g_viz_geodetic.m index f91a2a364017cfb8a8198ceb1b84b41e93fa3453..0b69b70be6ebf7c4b722cf0dd95890b1cf94ca68 100644 --- a/src/g_rect/g_viz_field.m +++ b/src/g_rect/g_viz_geodetic.m @@ -1,4 +1,4 @@ -function [h_figure,h_pcolor,h_datetext] = g_viz_field(imgFname,rectFile,varargin) +function [h_figure,h_pcolor,h_datetext] = g_viz_geodetic(imgFname,rectFile,varargin) %G_VIZ_FIELD Generates a map with a georectified image % G_VIZ_FIELD Generates a map with the georectified image in imgFname % Converts in grayscale if needed.