Commit 47ace827 authored by Daniel Bourgault's avatar Daniel Bourgault

The cdoe can now deal with GCPs with eleation. Plus the code now comes with m_map.

parent 13a4218a
......@@ -33,20 +33,18 @@ dH = 0.0;
dtheta = 20.0;
% Order of the polynomial correction (0, 1 or 2)
polyOrder = 1;
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
54 719 -70.54500 47.335
99 661 -70.505 47.375
452 641 -70.435 47.389
429 633 -70.418 47.408
816 644 -70.393 47.368
\ No newline at end of file
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
......@@ -45,18 +45,18 @@ precision = 'double';
% Ground Control Points (GCP).
gcpData = true;
2999 226 0.500 0.00
1694 220 1.500 0.00
528 677 2.290 0.50
289 1231 2.290 1.00
3819 686 0.000 0.50
4018 1266 0.000 1.00
2566 1235 0.890 0.99
3806 2118 0.255 1.56
4131 1580 0.000 1.23
521 2548 1.910 1.82
248 2328 2.060 1.71
14 1868 2.290 1.46
2999 226 0.500 0.00 0
1694 220 1.500 0.00 0
528 677 2.290 0.50 0
289 1231 2.290 1.00 0
3819 686 0.000 0.50 0
4018 1266 0.000 1.00 0
2566 1235 0.890 0.99 0
3806 2118 0.255 1.56 0
4131 1580 0.000 1.23 0
521 2548 1.910 1.82 0
248 2328 2.060 1.71 0
14 1868 2.290 1.46 0
......
......@@ -4,7 +4,7 @@ function errGeoFit = g_error_geofit(cv,imgWidth,imgHeight,xp,yp,ic,jc,...
hfovGuess,lambdaGuess,phiGuess,HGuess,thetaGuess,...
dhfov,dlambda,dphi,dH,dtheta,...
LON0,LAT0,...
i_gcp,j_gcp,lon_gcp,lat_gcp,...
i_gcp,j_gcp,lon_gcp0,lat_gcp0,h_gcp,...
theOrder,field);
......@@ -20,7 +20,6 @@ end
[LON LAT] = g_pix2ll(xp,yp,imgWidth,imgHeight,ic,jc,...
hfov,lambda,phi,theta,H,LON0,LAT0,field);
% Calculate the error between ground control points (GCPs)
% and image control points once georectified.
%
......@@ -35,6 +34,9 @@ ngcp = length(i_gcp);
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);
for k = 1:ngcp
% Calculate the distance (i.e. error) between GCP and rectificed ICPs.
......@@ -63,17 +65,3 @@ if abs(phi - phiGuess) > dphi; errGeoFit = inf; end
if abs(H - HGuess) > dH; errGeoFit = inf; end
if abs(theta - thetaGuess) > dtheta; errGeoFit = inf; end
%fprintf('\n')
%fprintf('Field of view (fov): %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('RMS error (m): %f\n',errGeoFit)
%fprintf('\n')
%figure(100)
%hold on;
%plot(errGeoFit);
......@@ -61,7 +61,6 @@ vfov = 2*atand(tand(hfov/2)/aspectRatio);
fx = (imgWidth/2)/tand(hfov/2);
fy = (imgHeight/2)/tand(vfov/2);
% Subtract the principal point
x_p = xp - x_c + (jc);
y_p = yp - y_c + (ic);
......
function [lon_gcp,lat_gcp] = g_proj_GCP(LON0,LAT0,H,lon_gcp0,lat_gcp0,h_gcp,field)
%G_PROJ_GCP
%
% This function projects the GCPs that have an elevation greater than 0
% onto the water level of zero elevation. This is a trick in order to be
% able to deal with GCPs that have elevation. The principle can be
% illustrated in 1D like this. Consider the position x of a GCP that has an
% elevation h. Then, the projection (xp) of this GCP onto the plane that has
% zero elevation and along a line that joins the camera position of elevation
% H and the GCP is given by:
%
% xp = x * (H - h)/H
%
% INPUTS:
%
% LON0: The longitiude of the camera
% LAT0: The latitude of the camera
% H: The elevation of the camera
% lon_gcp0: The longitudes of the GCPs
% lat_gcp0: The latitudes 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
%
n_gcp = length(h_gcp);
for i = 1:n_gcp
if field == true
% The calculation for the distance between the GCPs and the camera
% is done on an elliptical earth and this takes a long calculation
% with the command m_idist from the m_map package. Therefore, the
% calculation is done only on the GCPs with no elevation, hence
% the 'if' condition here.
if h_gcp(i) > 0
[range,a12,a21] = m_idist(LON0,LAT0,lon_gcp0(i),lat_gcp0(i));
proj_factor = H/(H-h_gcp(i));
new_range = range*proj_factor;
[lon_gcp(i),lat_gcp(i),a21] = m_fdist(LON0,LAT0,a12,new_range);
lon_gcp(i) = lon_gcp(i) - 360;
else
% If there's no elevation, then there's nothing to do.
lon_gcp(i) = lon_gcp0(i);
lat_gcp(i) = lat_gcp0(i);
end
else
% For lab 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;
range = sqrt(dx^2 + dy^2);
proj_factor = H/(H-h_gcp(i));
new_range = range*proj_factor;
beta = atan2(dy,dx);
lon_gcp(i) = new_range*cos(beta) + LON0;
lat_gcp(i) = new_range*sin(beta) + LAT0;
end
end
\ No newline at end of file
......@@ -37,12 +37,23 @@ function g_rect
%
% LAT0: Same as LON0 but for the latitude.
%
% lon_gcp: A vector containing the longitude of each ground
% lon0_gcp: A vector containing the longitude of each ground
% control points (GCP). For the lab case this is the
% cartesian coordinate in meters.
% cartesian coordinate in meters. These GCPs may have
% elevations.
%
% lat0_gcp: Same as lon_gcp for latitude.
%
% lon_gcp: A vector containing the longitude of each ground
% control points (GCP) projected onto the water level
% i.e. at 0 m of elevation. For the lab case this is
% the cartesian coordinate in meters.
%
% lat_gcp: Same as lon_gcp for latitude.
%
% h_gcp: The elevation in meter of the GCPs. The elevation is
% 0 m if taken at water level.
%
% i_gcp: The horizontal index of the image ground control
% points.
%
......@@ -53,11 +64,12 @@ function g_rect
%
% phi: Camera tilt angle [degree].
%
% lambda: Camera dip angle [degree].
%
% H: The camera altitude relative to the water [m].
%
% theta: View angle clockwise from North [degree].
%
%
% errGeoFit: The rms error of the georectified image after
% geometrical transformation [m].
%
......@@ -71,6 +83,7 @@ function g_rect
% and calculations can always be done in double
% precision.
%
% Delta: The effective resolution of the georectified image.
%
% Other m-files required: Works best with the m_map package for
% visualization.
......@@ -83,7 +96,7 @@ function g_rect
% February 2013
%
%
%%
% The minimization is repeated nMinimize times where each time a random
% combination of the initial guesses is chosen within the given
% uncertainties provided by the user. This is becasue the algorithm often
......@@ -91,7 +104,6 @@ function g_rect
% that the minimum found is a true minimum within the uncertainties provided.
nMinimize = 50;
%% Read the parameter file
% Count the number of header lines before the ground control points (GCPs)
% The GCPs start right after the variable gcpData is set to true.
......@@ -114,12 +126,28 @@ while gcpData == false
end
fclose(fid);
%% Import the GCP data (4 column) at the end of the parameter file
gcp = importdata(inputFname,' ',nHeaderLine)
%% 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_gcp = gcp.data(:,3);
lat_gcp = gcp.data(:,4);
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);
% Get the image size
......@@ -132,7 +160,6 @@ if precision == 'single'
imgHeight = single(imgHeight);
end
%% Display information
fprintf('\n')
fprintf(' INPUT PARAMETERS\n')
......@@ -169,9 +196,12 @@ colormap(gray);
hold on
for i = 1:ngcp
plot(i_gcp(i),j_gcp(i),'r.');
text(i_gcp(i),j_gcp(i),num2str(i),'color','r','horizontalalignment','right');
text(i_gcp(i),j_gcp(i),[num2str(i),'(',num2str(h_gcp(i)),')'],...
'color','r',...
'horizontalalignment','right',...
'fontsize',16);
end
title('Ground Control Points','color','r');
title('Ground Control Points Number (elevation in meters)','color','r');
print -dpng -r300 image1.png
......@@ -179,11 +209,9 @@ fprintf('\n')
ok = input('Is it ok to proceed with the rectification (y/n): ','s');
if ok ~= 'y'
return
%break
end
%%
nUnknown = 0;
if dhfov > 0.0; nUnknown = nUnknown+1; end
if dlambda > 0.0; nUnknown = nUnknown+1; end
......@@ -196,7 +224,6 @@ if nUnknown > ngcp
fprintf('WARNING: \n');
fprintf(' The number of GCPs is < number of unknown parameters.\n');
fprintf(' Program stopped.\n');
%break
return
end
......@@ -227,7 +254,6 @@ if nUnknown > 0
%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;
......@@ -297,10 +323,9 @@ if nUnknown > 0
hfovGuess,lambdaGuess,phiGuess,HGuess,thetaGuess,...
dhfov,dlambda,dphi,dH,dtheta,...
LON0,LAT0,...
i_gcp,j_gcp,lon_gcp,lat_gcp,...
i_gcp,j_gcp,lon_gcp0,lat_gcp0,h_gcp,...
theOrder,field);
if errGeoFit < bestErrGeoFit
bestErrGeoFit = errGeoFit;
cvBest = cv;
......@@ -334,9 +359,11 @@ if nUnknown > 0
if length(theOrder) > 1
fprintf('The rms error after geometrical correction (m): %f\n',bestErrGeoFit);
end
end
% Project the GCP that have elevation.
[lon_gcp,lat_gcp] = g_proj_GCP(LON0,LAT0,H,lon_gcp0,lat_gcp0,h_gcp,field);
%%
% Now construct the matrices LON and LAT for the entire image using the
......@@ -369,7 +396,8 @@ fprintf('Saving rectification file in: %s\n',outputFname);
save(outputFname,'imgFname','firstImgFname','lastImgFname',...
'LON','LAT',...
'LON0','LAT0',...
'lon_gcp','lat_gcp',...
'lon_gcp0','lat_gcp0',...
'lon_gcp','lat_gcp','h_gcp',...
'i_gcp','j_gcp',...
'hfov','lambda','phi','H','theta',...
'errGeoFit','errPolyFit',...
......
......@@ -18,6 +18,12 @@ function [h_figure,h_pcolor,h_datetext] = g_viz_field(imgFname,rectFile,varargin
% corresponding objects in the figure, for use in g_viz_anim
%
% Set some plotting parameters
ms = 10; % Marker Size
fs = 10; % Font size
lw = 2; % Line width
% Option
show_time = 0;
show_land = 0;
land_color = [241 235 144]/255;
......@@ -34,30 +40,21 @@ if length(varargin) > 1
end
end
% Set some plotting parameters
ms = 10; % Marker Size
fs = 10; % Font size
lw = 2; % Line width
% Load the georectification file
load(rectFile);
%p = size(LON,3);
lon_min = min(lon_gcp);
lon_max = max(lon_gcp);
lat_min = min(lat_gcp);
lat_max = max(lat_gcp);
lon_min = min(lon_min,LON0);
lon_max = max(lon_max,LON0);
lat_min = min(lat_min,LAT0);
lat_max = max(lat_max,LAT0);
% Determine the region to plot, delimited by the GCP and the camera
% position. Add a factor (fac) all around
fac = 0.1;
lon_min= lon_min - fac*abs(lon_max-lon_min);
lon_max= lon_max + fac*abs(lon_max-lon_min);
lat_min= lat_min - fac*abs(lat_max-lat_min);
lat_max= lat_max + fac*abs(lat_max-lat_min);
lon_min = min([lon_gcp,LON0]);
lon_max = max([lon_gcp,LON0]);
lat_min = min([lat_gcp,LAT0]);
lat_max = max([lat_gcp,LAT0]);
lon_min = lon_min - fac*abs(lon_max - lon_min);
lon_max = lon_max + fac*abs(lon_max - lon_min);
lat_min = lat_min - fac*abs(lat_max - lat_min);
lat_max = lat_max + fac*abs(lat_max - lat_min);
rgb0 = double(imread(imgFname))/255;
......@@ -70,19 +67,19 @@ end
clear rgb0;
int = int';
int = int - nanmean(nanmean(int));
% Normalize the image if wanted
%int = int - nanmean(nanmean(int));
figure('Renderer', 'painters', 'Position', [100 100 900 700])
%figure('Renderer', 'painters', 'Position', [100 100 900 700])
figure
m_proj('mercator','longitudes',[lon_min lon_max],'latitudes',[lat_min lat_max]);
hold on;
[X,Y] = m_ll2xy(LON,LAT);
cmap = contrast(int,256);
colormap(cmap);
h = pcolor(X,Y,int);
h = m_pcolor(LON,LAT,int);
shading('interp');
% Uncomment one of these lines if you want the coastline to be plotted.
if show_land
if strcmpi(show_land,'fjord')
......@@ -109,7 +106,14 @@ m_plot(LON0,LAT0,'kx','markersize',ms,'linewidth',lw); % Camera location
%% Plot GCPs and ICPs.
for n=1:length(i_gcp)
% Plot the original GCPs that may have elevations
m_plot(lon_gcp0(n),lat_gcp0(n),'ko','markersize',ms,'linewidth',lw);
% Plot the projected GCPs that are actually used for the minimization
m_plot(lon_gcp(n),lat_gcp(n),'bo','markersize',ms,'linewidth',lw);
% Plot the Image Control Points once georectified
m_plot(LON(i_gcp(n),j_gcp(n)),LAT(i_gcp(n),j_gcp(n)),'rx','MarkerSize',ms,'linewidth',lw);
end
......@@ -118,7 +122,7 @@ end
clear LON LAT X Y
m_grid_old('box','fancy','fontsize',fs);
m_grid('box','fancy','fontsize',fs);
%m_gshhs('f','patch', 'gray');
%print('g_rect_drone_1', '-dpng')
if nargout > 1
......
function g_viz_lab(imgFname,outputFname);
function g_viz_lab(imgFname,rectFile);
% Set some plotting parameters
ms = 10; % Marker Size
fs = 10; % Font size
lw = 2; % Line width
load(outputFname);
lon_min = min(lon_gcp);
lon_max = max(lon_gcp);
lat_min = min(lat_gcp);
lat_max = max(lat_gcp);
lon_min = min(lon_min,LON0);
lon_max = max(lon_max,LON0);
lat_min = min(lat_min,LAT0);
lat_max = max(lat_max,LAT0);
load(rectFile);
% Determine the region to plot, delimited by the GCP and the camera
% position. Add a factor (fac) all around
fac = 0.1;
lon_min= lon_min - fac*abs(lon_max-lon_min);
lon_max= lon_max + fac*abs(lon_max-lon_min);
lat_min= lat_min - fac*abs(lat_max-lat_min);
lat_max= lat_max + fac*abs(lat_max-lat_min);
lon_min = min([lon_gcp,LON0]);
lon_max = max([lon_gcp,LON0]);
lat_min = min([lat_gcp,LAT0]);
lat_max = max([lat_gcp,LAT0]);
lon_min = lon_min - fac*abs(lon_max - lon_min);
lon_max = lon_max + fac*abs(lon_max - lon_min);
lat_min = lat_min - fac*abs(lat_max - lat_min);
lat_max = lat_max + fac*abs(lat_max - lat_min);
rgb0 = double(imread(imgFname))/255;
......@@ -35,18 +33,23 @@ int = int';
hold on;
X = LON;
Y = LAT;
colormap(gray);
h = pcolor(X,Y,int);
pcolor(LON,LAT,int);
shading('flat');
plot(LON0,LAT0,'kx','markersize',ms,'linewidth',lw); % Camera location
%% Plot GCPs and ICPs.
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);
% Plot the projected GCPs that are actually used for the minimization
plot(lon_gcp(n),lat_gcp(n),'bo','markersize',ms,'linewidth',lw);
% 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
axis([lon_min lon_max lat_min lat_max]);
\ No newline at end of file
function g_roi2(imgFname)
% G_ROI is a function to determine regions of interest in an image.
%
% This function allows the user to define a series of polygons
% on an image in order to define regions of interest (roi) for
% image stabilization. The result is a 2D matrix roi with 1s and 0s.
% Just follow the instructions.
%
% Input: imgFname: the image filename
% Output: The 2D matrix roi written to filename roi.mat
%
% ** There is tips to modifie and adjuste the roi before saving it in the
% help menu of the impoly function.
% Example : Hold the button A and click on the polygon border to add
% points.
%
% Author: Daniel Bourgault - 2011
%
% Updated : - 11/05/19 - Show the lines while creating the polygon.
% Updated : - 08/08/19 - Enable modification of the polygon before saving.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Reading reference images
im = imread(imgFname);
[m,n,~] = size(im);
roi(1:m,1:n) = 0; % Create ROI matrix
%% Generating Region of Interest
imagesc(im)
hold on
colormap(gray);
continue_draw = 'y';
while continue_draw == 'y'
% Instructions
fprintf(' Make a polygon with left click on the mouse\n')
fprintf(' Right click to link the last and first point\n')
fprintf(' Double click on the polygon when you are done arranging points\n')
fprintf(' You will have the option to make more polygons when done with this one\n')
% Drawing
h=impoly;
wait(h);
Mask = createMask(h);
roi = roi+Mask;
continue_draw = input('Draw another polygon (y/n) ? ','s');
end
hold off
ij = roi > 1;
roi(ij) = 1;
% Final figure
figure(2)
imagesc(roi)
save([imgFname(1:end-4),'_roi.mat'],'roi','imgFname');
\ No newline at end of file
% M_Map - mapping toolbox (Author: rich@eos.ubc.ca)
% Version 1.4m Feb 2020
%
% You have collected your data, loaded it into Matlab, analyzed
% everything to death, and now you want to make a simple map showing
% how it relates to the world.
%
% But you can't.
%