g_error_geofit.m 2.46 KB
Newer Older
Daniel Bourgault's avatar
Daniel Bourgault committed
1
2
3
4
5
6
function errGeoFit = g_error_geofit(cv,imgWidth,imgHeight,xp,yp,ic,jc,...
                                   hfov,lambda,phi,H,theta,...
                                   hfov0,lambda0,phi0,H0,theta0,...
                                   hfovGuess,lambdaGuess,phiGuess,HGuess,thetaGuess,...
                                   dhfov,dlambda,dphi,dH,dtheta,...
                                   LON0,LAT0,...
7
                                   i_gcp,j_gcp,lon_gcp0,lat_gcp0,h_gcp,...
Daniel Bourgault's avatar
Daniel Bourgault committed
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
                                   theOrder,field);

                               
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
  if theOrder(i) == 4; H      = cv(i); end
  if theOrder(i) == 5; theta  = cv(i); end
end

% Perform the geometrical transformation
[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. 
%
% This error (errGeoFit) is the error associated with the geometrical fit,
% as opposed to the error associated with a polynomial fit that can be 
% calculated if requested by the user. 
%

% Determine the number of CGP
ngcp = length(i_gcp);

ngcpFinite = 0;
errGeoFit  = 0.0;

37
38
% 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);
Daniel Bourgault's avatar
Daniel Bourgault committed
39

40
41
for k = 1:ngcp
    
Daniel Bourgault's avatar
Daniel Bourgault committed
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
  % 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);
  
  % Check if the distance is finite. The distance may be NaN for some
  % GCPs that may temporarily be above the horizon. Those points are 
  % blanked out in the function g_pix2ll.
  if isfinite(distance) == 1 
    errGeoFit = errGeoFit + distance^2;
    ngcpFinite = ngcpFinite + 1;
  else
    errGeoFit = Inf;
  end
  
end

% rms distance
errGeoFit = sqrt(errGeoFit/ngcpFinite);

% Check if the parameters are within the specified uncertainties.
% If not set the error to infinity.
if abs(hfov - hfovGuess)     > dhfov;   errGeoFit = inf; end
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