i.ortho.photo (#401)

Add optional correction for panorama cameras, adjust for earth curvature

* i.ortho.elev

improve error messages, set mapset from target if not given

* i.ortho.photo

Omega and Phi are swapped in the description but correct in the code

* i.ortho

support panorama corrections

* i.ortho photo2image

InsertColumnInfo -> InsertColumn

* i.ortho

support orthorectification without GCPs if camera position and angles are known

* i.ortho.transform

fix for orthorectification

* i.ortho.photo

fix g.gui.image2target for wxpython 4 (phoenix)

* i.ortho.target

XY and ll locations are not supported

* i.ortho.init

explain importance of standard deviations for the improvement of the orthorectification parameters

* i.ortho.photo

lib: report number of iterations needed for the refinement of orthorectification parameters

* i.ortho.photo

refine panorama correction
add correction for earth curvature

* i.ortho.photo

I_compute_ortho_equations() must be run in the target location for earth curvature adjustments

gui/wxpython/image2target/ii2t_manager.py

@@ -2241,7 +2241,7 @@ class GCPList(ListCtrl,
                         _('Forward error'),
                         _('Backward error')):
                 info.SetText(lbl)
-                self.InsertColumnInfo(idx_col, info)
+                self.InsertColumn(idx_col, info)
                 idx_col += 1
 
     def LoadData(self):

gui/wxpython/photo2image/ip2i_manager.py

@@ -1554,7 +1554,7 @@ class GCPList(ListCtrl,
                         _('Forward error'),
                         _('Backward error')):
                 info.SetText(lbl)
-                self.InsertColumnInfo(idx_col, info)
+                self.InsertColumn(idx_col, info)
                 idx_col += 1
 
     def LoadData(self):

imagery/i.ortho.photo/i.ortho.elev/main.c

@@ -48,13 +48,13 @@ int main(int argc, char *argv[])
     char group[GNAME_MAX];
 
     char *elev_layer;
-    char *mapset_elev;
-    char *location_elev;
+    char *mapset_elev, *mapset_elev_old;
+    char *location_elev, *location_elev_old;
     char *math_exp;
     char *units;
     char *nd;
 
-    char buf[100];
+    char buf[GPATH_MAX];
     int stat;
     int overwrite;
 
@@ -111,30 +111,34 @@ int main(int argc, char *argv[])
 
     elev_layer = (char *)G_malloc(GNAME_MAX * sizeof(char));
     mapset_elev = (char *)G_malloc(GMAPSET_MAX * sizeof(char));
+    mapset_elev_old = (char *)G_malloc(GMAPSET_MAX * sizeof(char));
     location_elev = (char *)G_malloc(80 * sizeof(char));
+    location_elev_old = (char *)G_malloc(80 * sizeof(char));
     math_exp = (char *)G_malloc(80 * sizeof(char));
     units = (char *)G_malloc(80 * sizeof(char));
     nd = (char *)G_malloc(80 * sizeof(char));
 
     *elev_layer = 0;
     *mapset_elev = 0;
+    *mapset_elev_old = 0;
     *location_elev = 0;
+    *location_elev_old = 0;
     *math_exp = 0;
     *units = 0;
     *nd = 0;
 
     strcpy(group, group_opt->answer);
-    if(loc_opt->answer)
+    if (loc_opt->answer)
     	strcpy(location_elev, loc_opt->answer);
-    if(mapset_opt->answer)
+    if (mapset_opt->answer)
     	strcpy(mapset_elev, mapset_opt->answer);
     /*if(elev_opt->answer)
     	strcpy(elev_layer, elev_opt->answer);*/
-    if(math_opt->answer)
+    if (math_opt->answer)
     	strcpy(math_exp, math_opt->answer);
-    if(unit_opt->answer)
+    if (unit_opt->answer)
     	strcpy(units, unit_opt->answer);
-    if(nd_opt->answer)
+    if (nd_opt->answer)
     	strcpy(nd, nd_opt->answer);
 	
     if (!I_get_target(group, location, mapset)) {
@@ -149,17 +153,18 @@ int main(int argc, char *argv[])
     /*Report the contents of the ELEVATION file as in the GROUP */
     if (print_flag->answer) {
 	/*If the content is empty report an error */
-	if(!I_get_group_elev(group, elev_layer, mapset_elev, location_elev, math_exp, units, nd)){
+	if (!I_get_group_elev(group, elev_layer, mapset_elev, location_elev, math_exp, units, nd)) {
 		G_fatal_error(_("Cannot find default elevation map for target in group [%s]"),group);
-	/*If there is a content, report it as a message */
-	} else {
-		G_message("map:\t\t\t%s",elev_layer);
-		G_message("mapset:\t\t\t%s",mapset_elev);
-		G_message("location:\t\t%s",location_elev);
-		G_message("math expression:\t%s",math_exp);
-		G_message("units:\t\t\t%s",units);
-		G_message("nodata value:\t\t%s",nd);
-		exit(EXIT_SUCCESS);
+	}
+	/*If there is a content, print it */
+	else {
+	    fprintf(stdout, "map:\t\t\t%s\n",elev_layer);
+	    fprintf(stdout, "mapset:\t\t\t%s\n",mapset_elev);
+	    fprintf(stdout, "location:\t\t%s\n",location_elev);
+	    fprintf(stdout, "math expression:\t%s\n",math_exp);
+	    fprintf(stdout, "units:\t\t\t%s\n",units);
+	    fprintf(stdout, "nodata value:\t\t%s\n",nd);
+	    exit(EXIT_SUCCESS);
 	}    
     }
 
@@ -190,27 +195,44 @@ int main(int argc, char *argv[])
 	    G_fatal_error(_("Elevation map name is missing. Please set '%s' option"),
 	        elev_opt->key);
 	}
-	/* elevation map exists in target ? */
-	if (G_find_raster2(elev_opt->answer, mapset) == NULL) {
-            /* select current location */
-	    select_current_env();
-	    G_fatal_error(_("Raster map <%s> not found"), elev_opt->answer);
-	}
 	
 	/* return to current Location/mapset to write in the group file */
 	select_current_env();
 	
 	/* load information from the ELEVATION file in the GROUP */
-	I_get_group_elev(group, elev_layer, mapset_elev, location_elev, math_exp, units, nd);
-	/* Modify ELEVATION file in source GROUP */
-	I_put_group_elev(group,elev_opt->answer,mapset_opt->answer,loc_opt->answer, 
-			math_opt->answer, unit_opt->answer, nd_opt->answer);
+	if (I_find_group_elev_file(group)) {
+	    I_get_group_elev(group, elev_layer, mapset_elev_old, location_elev_old,
+	                     math_exp, units, nd);
+	    if (*location_elev == 0)
+		strcpy(location_elev, location_elev_old);
+	    if (*mapset_elev == 0)
+		strcpy(mapset_elev, mapset_elev_old);
+	}
+	/* if location and/or mapset of elevation are not set, use target */
+	if (*mapset_elev == 0)
+	    strcpy(mapset_elev, mapset);
+	if (*location_elev == 0)
+	    strcpy(location_elev, location);
+
+	/* select target location */
+	select_target_env();
+	/* elevation map exists in target ? */
+	if (G_find_raster2(elev_opt->answer, mapset_elev) == NULL) {
+            /* select current location */
+	    select_current_env();
+	    G_fatal_error(_("Raster map <%s> not found"), elev_opt->answer);
+	}
 	/* select current location */
 	select_current_env();
 
+	/* Modify ELEVATION file in source GROUP */
+	I_put_group_elev(group, elev_opt->answer, mapset_elev, location_elev, 
+			 math_exp, units, nd);
+
         G_message(_("Group [%s] in location [%s] mapset [%s] now uses elevation map [%s]"),
-	          group, location, mapset, elev_opt->answer);
-    }else{
+	          group, G_location(), G_mapset(), elev_opt->answer);
+    }
+    else {
 	G_fatal_error(_("Mapset [%s] in target location [%s] - %s "),
                   mapset, location,
 		  stat == 0 ? _("permission denied\n") : _("not found\n"));

imagery/i.ortho.photo/i.ortho.init/i.ortho.init.html

@@ -32,15 +32,15 @@ The following menu is displayed:
 	INITIAL XC: Meters                __________
 	INITIAL YC: Meters                __________
 	INITIAL ZC: Meters                __________
-	INITIAL omega (roll) degrees:     __________
-	INITIAL phi  (pitch) degrees:     __________
+	INITIAL omega (pitch) degrees:    __________
+	INITIAL phi  (roll) degrees:      __________
 	INITIAL kappa  (yaw) degrees:     __________
 
 	Standard Deviation XC: Meters     __________
 	Standard Deviation YC: Meters     __________
 	Standard Deviation ZC: Meters     __________
-	Std. Dev. omega (roll) degrees:   __________
-	Std. Dev. phi  (pitch) degrees:   __________
+	Std. Dev. omega (pitch) degrees:  __________
+	Std. Dev. phi  (roll) degrees:    __________
 	Std. Dev. kappa  (yaw) degrees:   __________
 
         Use these values at run time? (1=yes, 0=no)
@@ -67,24 +67,33 @@ represent an approximation for the cameras attitude  at the time of
 exposure. 
 
 <ul>
-<li> Omega (roll): Raising or lowering of the wings (turning around the
-  aircraft's axis);
-<li> Phi (pitch): Raising or lowering of the aircraft's front (turning
+<li> Omega (pitch): Raising or lowering of the aircraft's front (turning
   around the wings' axis);
+<li> Phi (roll): Raising or lowering of the wings (turning around the
+  aircraft's axis);
 <li> Kappa (yaw): Rotation needed to align the aerial photo to true north:
   needs to be denoted as +90 degree for clockwise turn and -90 degree for
   a counterclockwise turn.
 </ul>
+
+<p>
+If ground control points are available, the INITIAL values are iteratively 
+corrected. This is particularl useful when the INITIAL values are rather 
+rough estimates.
+
 <p>
 
 The standard deviations for (XC,YC,ZC) are expressed in meters, and
 are used as <em>a priori</em> values for the standard deviations used in 
-computation of the ortho rectification parameters.
+computation of the ortho rectification parameters. Higher values improve 
+the refinement of the initial camera exposure. As a rule of thumb, 5% 
+of the estimated target extents should be used.
 <p>
 
 The standard deviations for (omega,phi,kappa) are expressed in degrees, and
 are used as <em>a priori</em> values for the standard deviations used in 
-computation of the ortho rectification parameters.
+computation of the ortho rectification parameters. As a rule of thumb, 
+2 degrees should be used.
 
 <p>
 If <i>Use these values at run time? (1=yes, 0=no)</i> is set to 0, the

imagery/i.ortho.photo/i.ortho.init/main.c

@@ -89,12 +89,12 @@ int main(int argc, char *argv[])
     omega_opt = G_define_option();
     omega_opt->key = "omega";
     omega_opt->type = TYPE_DOUBLE;
-    omega_opt->description = _("Initial Camera Omega (roll) degrees");
+    omega_opt->description = _("Initial Camera Omega (pitch) degrees");
 
     phi_opt = G_define_option();
     phi_opt->key = "phi";
     phi_opt->type = TYPE_DOUBLE;
-    phi_opt->description = _("Initial Camera Phi (pitch) degrees");
+    phi_opt->description = _("Initial Camera Phi (roll) degrees");
 
     kappa_opt = G_define_option();
     kappa_opt->key = "kappa";
@@ -104,12 +104,12 @@ int main(int argc, char *argv[])
     omegasd_opt = G_define_option();
     omegasd_opt->key = "omega_sd";
     omegasd_opt->type = TYPE_DOUBLE;
-    omegasd_opt->description = _("Apriori Omega (roll) standard deviation");
+    omegasd_opt->description = _("Apriori Omega (pitch) standard deviation");
 
     phisd_opt = G_define_option();
     phisd_opt->key = "phi_sd";
     phisd_opt->type = TYPE_DOUBLE;
-    phisd_opt->description = _("Apriori Phi (pitch) standard deviation");
+    phisd_opt->description = _("Apriori Phi (roll) standard deviation");
 
     kappasd_opt = G_define_option();
     kappasd_opt->key = "kappa_sd";

imagery/i.ortho.photo/i.ortho.photo/i.ortho.photo.html

@@ -227,13 +227,13 @@ may be selected or modified.
     <li><b>X</b>: East aircraft position;</li>
     <li><b>Y</b>: North aircraft position;</li>
     <li><b>Z</b>: Flight heigh above surface;</li>
-	<li><b>Omega (roll)</b>: Raising or lowering of the wings (turning 
-	around the aircraft's axis);</li>
-    <li><b>Phi (pitch)</b>: Raising or lowering of the aircraft's front 
+    <li><b>Omega (pitch)</b>: Raising or lowering of the aircraft's front 
 	(turning around the wings' axis);</li>
-	<li><b>Kappa (yaw)</b>: Rotation needed to align the aerial photo to 
-	true north: needs to be denoted as +90&deg; for clockwise turn and
-	-90&deg; for a counter-clockwise turn.</li>
+    <li><b>Phi (roll)</b>: Raising or lowering of the wings (turning 
+    around the aircraft's axis);</li>
+    <li><b>Kappa (yaw)</b>: Rotation needed to align the aerial photo to 
+    true north: needs to be denoted as +90&deg; for clockwise turn and
+    -90&deg; for a counter-clockwise turn.</li>
 </ul>
 
 <div align="center" style="margin: 10px">
@@ -325,6 +325,12 @@ and further image analysis.
 </li>
 </ol>
 
+<h2>REFERENCES</h2>
+Wolf P.R. (1983). <i>Elements of Photogrammetry: With Air Photo 
+Interpretation and Remote Sensing</i>
+<b>McGraw Hill Higher Education</b>
+ISBN-10: 0070713456, ISBN-13: 978-0070713451 <br>
+
 <h2>SEE ALSO</h2>
 
 <em>

imagery/i.ortho.photo/i.ortho.rectify/cp.c

@@ -4,17 +4,19 @@
 
 int get_conz_points(struct Ortho_Image_Group *group)
 {
-    char msg[200];
+    char msg[1024];
 
     if (!I_get_con_points(group->name, &group->control_points))
-	exit(0);
+	group->control_points.count = 0;
 
     sprintf(msg, _("Control Z Point file for group <%s@%s> - "),
 	    group->name, G_mapset());
 
     G_verbose_message(_("Computing equations..."));
 
+    select_target_env();
     Compute_ortho_equation(group);
+    select_current_env();
 
     switch (group->con_equation_stat) {
     case -1:
@@ -37,7 +39,7 @@ int get_conz_points(struct Ortho_Image_Group *group)
 
 int get_ref_points(struct Ortho_Image_Group *group)
 {
-    char msg[200];
+    char msg[1024];
 
     if (!I_get_ref_points(group->name, &group->photo_points))
 	exit(0);

imagery/i.ortho.photo/i.ortho.rectify/main.c

@@ -78,7 +78,7 @@ int main(int argc, char *argv[])
      *interpol;			/* interpolation method:
 				   nearest neighbor, bilinear, cubic */
 
-    struct Flag *c, *a;
+    struct Flag *c, *a, *pan_flag;
     struct GModule *module;
 
     G_gisinit(argv[0]);
@@ -139,6 +139,10 @@ int main(int argc, char *argv[])
     a->key = 'a';
     a->description = _("Rectify all raster maps in group");
 
+    pan_flag = G_define_flag();
+    pan_flag->key = 'p';
+    pan_flag->description = _("Enable panorama camera correction");
+
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
@@ -256,15 +260,19 @@ int main(int argc, char *argv[])
 		      group.name);
     }
 
+    /* panorama camera correction */
+    if (pan_flag->answer)
+	I_ortho_panorama();
+
+    /* get the target */
+    get_target(group.name);
+
     /* read the reference points for the group, compute image-to-photo trans. */
     get_ref_points(&group);
 
     /* read the control points for the group, convert to photo coords. */
     get_conz_points(&group);
 
-    /* get the target */
-    get_target(group.name);
-
     /* Check the GRASS_OVERWRITE environment variable */
     if ((overstr = getenv("GRASS_OVERWRITE")))  /* OK ? */
 	target_overwrite = atoi(overstr);

imagery/i.ortho.photo/i.ortho.target/main.c

@@ -40,6 +40,8 @@ int main(int argc, char *argv[])
     /* newly defined location and maspet */
     char target_location[GMAPSET_MAX];
     char target_mapset[GMAPSET_MAX];
+    int stat;
+    struct Cell_head target_window;
 
     G_gisinit(argv[0]);
 
@@ -75,6 +77,20 @@ int main(int argc, char *argv[])
 
     I_get_target(group, location, mapset);
     G_create_alt_env();
+    G_setenv_nogisrc("LOCATION_NAME", target_location);
+    stat = G_mapset_permissions(target_mapset);
+    if (stat != 1) {
+	G_fatal_error(_("Unable to access target location/mapset %s/%s"),
+	              target_location, target_mapset);
+    }
+
+    G_setenv_nogisrc("MAPSET", target_mapset);
+    G_get_window(&target_window);
+    if (target_window.proj == PROJECTION_XY)
+	G_fatal_error(_("Target locations with XY (unreferenced) are not supported"));
+    else if (target_window.proj == PROJECTION_LL)
+	G_fatal_error(_("Target locations with lon/lat are not supported"));
+
     G_switch_env();
     I_put_target(group, target_location, target_mapset);
 

imagery/i.ortho.photo/i.ortho.transform/main.c

@@ -16,6 +16,7 @@
  *
  *****************************************************************************/
 
+#include <unistd.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
@@ -55,6 +56,72 @@ static struct Ortho_Control_Points *points;
 static int count;
 static struct Stats fwd, rev;
 
+/* target fns taken from i.ortho.rectify */
+static int which_env = -1;	/* 0 = cur, 1 = target */
+
+int select_current_env(void)
+{
+    if (which_env < 0) {
+	G_create_alt_env();
+	which_env = 0;
+    }
+    if (which_env != 0) {
+	G_switch_env();
+	which_env = 0;
+    }
+
+    return 0;
+}
+
+int select_target_env(void)
+{
+    if (which_env < 0) {
+	G_create_alt_env();
+	which_env = 1;
+    }
+    if (which_env != 1) {
+	G_switch_env();
+	which_env = 1;
+    }
+
+    return 0;
+}
+
+static int get_target(void)
+{
+    char location[GMAPSET_MAX];
+    char mapset[GMAPSET_MAX];
+    char buf[1024];
+    int stat;
+
+    if (!I_get_target(group.name, location, mapset)) {
+	sprintf(buf, _("Target information for group <%s> missing"), group.name);
+	goto error;
+    }
+
+    sprintf(buf, "%s/%s", G_gisdbase(), location);
+    if (access(buf, 0) != 0) {
+	sprintf(buf, _("Target location <%s> not found"), location);
+	goto error;
+    }
+    select_target_env();
+    G_setenv_nogisrc("LOCATION_NAME", location);
+    stat = G_mapset_permissions(mapset);
+    if (stat > 0) {
+	G_setenv_nogisrc("MAPSET", mapset);
+	select_current_env();
+	return 1;
+    }
+    sprintf(buf, _("Mapset <%s> in target location <%s> - "), mapset, location);
+    strcat(buf, stat == 0 ? _("permission denied") : _("not found"));
+  error:
+    strcat(buf, "\n");
+    strcat(buf, _("Please run i.target for group "));
+    strcat(buf, group.name);
+    G_fatal_error("%s", buf);
+    return 1;			/* never reached */
+}
+
 static void update_max(struct Max *m, int n, double k)
 {
     if (k > m->val) {
@@ -375,13 +442,13 @@ static void do_pt_xforms(void)
     	fclose(fp);
 }
 
-
 int main(int argc, char **argv)
 {
     struct Option *grp, *fmt, *xfm_pts;
-    struct Flag *sum, *rev_flag, *dump_flag;
+    struct Flag *sum, *rev_flag, *dump_flag, *pan_flag;
     struct GModule *module;
     char *desc;
+    char *camera;
 
     G_gisinit(argv[0]);
 
@@ -438,6 +505,10 @@ int main(int argc, char **argv)
     dump_flag->key = 'x';
     dump_flag->description = _("Display transform matrix coefficients");
 
+    pan_flag = G_define_flag();
+    pan_flag->key = 'p';
+    pan_flag->description = _("Enable panorama camera correction");
+
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
@@ -450,20 +521,44 @@ int main(int argc, char **argv)
     forward = !rev_flag->answer;
     coord_file = xfm_pts->answer;
 
+    if (pan_flag->answer)
+	I_ortho_panorama();
+
     if (!I_get_ref_points(group.name, &group.photo_points)) {
 	G_fatal_error(_("Can not read reference points for group <%s>"),
 	              group.name);
     }
     if (!I_get_con_points(group.name, &group.control_points)) {
-	G_fatal_error(_("Can not read control points for group <%s>"),
-	              group.name);
+	group.control_points.count = 0;
     }
     
+    camera = (char *)G_malloc(GNAME_MAX * sizeof(char));
+    if (!I_get_group_camera(group.name, camera))
+	G_fatal_error(_("No camera reference file selected for group <%s>"),
+		      group.name);
+
+    if (!I_get_cam_info(camera, &group.camera_ref))
+	G_fatal_error(_("Bad format in camera file for group <%s>"),
+		      group.name);
+
+    /* get initial camera exposure station, if any */
+    if (I_find_initial(group.name)) {
+	if (!I_get_init_info(group.name, &group.camera_exp))
+	    G_warning(_("Bad format in initial exposure station file for group <%s>"),
+		      group.name);
+    }
+
+    /* get the target */
+    get_target();
+    
     points = &group.control_points;
 
     parse_format();
 
+    /* I_compute_ortho_equations() must be run in the target location */
+    select_target_env();
     compute_transformation();
+    select_current_env();
 
     analyze();
 

imagery/i.ortho.photo/lib/orthophoto.h

@@ -141,6 +141,8 @@ int I_ortho_ref(double, double, double, double *, double *, double *,
 int I_inverse_ortho_ref(double, double, double, double *, double *, double *,
 			struct Ortho_Camera_File_Ref *, double, double,
 			double, MATRIX);
+void I_ortho_panorama(void);
+
 /* ref_points.c */
 int I_new_ref_point(struct Ortho_Photo_Points *, double,
 		    double, double, double, int);

imagery/i.ortho.photo/lib/orthoref.c

@@ -44,6 +44,15 @@ FILE *debug;
 char msg[120];
 #endif
 
+static int panorama = 0;
+static double ellps_a = 0;
+
+/* enable panorama camera correction, 
+ * e.g. for CORONA KH-4A/B */
+void I_ortho_panorama(void)
+{
+    panorama = 1;
+}
 
 /* Compute the ortho rectification parameters */
 /* Camera position: XC, YC, ZC */
@@ -65,10 +74,15 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
     double Q1;
     double kappa1, kappa2, XC_var, YC_var, ZC_var;
     double omega_var, phi_var, kappa_var;
-    int i, iter, n;
+    int i, iter, max_iters, n;
     int first, active;
+    double e_a, e2;
 
     Q1 = (double)1.0;
+    
+    if (!G_get_ellipsoid_parameters(&e_a, &e2))
+	G_fatal_error(_("No ellpsoid parameters available"));
+    ellps_a = e_a;
 
     /* DEBUG */
 #ifdef DEBUG
@@ -79,17 +93,14 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
     }
 #endif
 
-    /*  Need at least 4 active control points */
+    /*  Need at least 4 active control points, 
+     *  or user-provided camera position */
     active = 0;
-    for (i = 0; i < cpz->count; i++) {
-	if (cpz->status[i] > 0)
-	    active++;
-    }
-    if (active < 4) {
-#ifdef DEBUG
-	fclose(debug);
-#endif
-	return (-1);
+    if (cpz) {
+	for (i = 0; i < cpz->count; i++) {
+	    if (cpz->status[i] > 0)
+		active++;
+	}
     }
 
     /*  Initialize (and zero out) all matrices needed */
@@ -169,6 +180,7 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
      * - init_info->status allows for disactivation of defined values 
      */
     if ((init_info != NULL) && (init_info->status == 1)) {
+	G_verbose_message("Using provided values for XC,YC,ZC,omega,phi,kappa");
 	/* Have initial values */
 	*XC = init_info->XC_init;
 	*YC = init_info->YC_init;
@@ -185,10 +197,11 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
 	WT1.x[4][4] = (Q1 / (init_info->phi_var * init_info->phi_var));
 	WT1.x[5][5] = (Q1 / (init_info->kappa_var * init_info->kappa_var));
     }
-    else {  /* set from mean values of active control points */
+    else if (active >= 4) {  /* set from mean values of active control points */
 	double dist_grnd, dist_photo;
 	double meanx, meany, meanX, meanY, meanZ;
 	
+	G_verbose_message("Estimating values for XC,YC,ZC,omega,phi,kappa");
 	/* set initial XC and YC from mean values of control points */
 	meanx = meany = meanX = meanY = meanZ = 0;
 	x = y = X = Y = 0;
@@ -263,9 +276,18 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
 	WT1.x[4][4] = (Q1 / (phi_var * phi_var));
 	WT1.x[5][5] = (Q1 / (kappa_var * kappa_var));
     }
-    G_debug(1, "INITIAL CAMERA POSITION:");
-    G_debug(1, "XC: %.2f, YC: %.2f, ZC: %.2f", *XC, *YC, *ZC);
-    G_debug(1, "Omega %.2f, Phi %.2f, Kappa: %.2f", *Omega, *Phi, *Kappa);
+    else {
+	G_warning(_("Camera postion not available"));
+#ifdef DEBUG
+	fclose(debug);
+#endif
+	return (0);
+    }
+
+    G_verbose_message("INITIAL CAMERA POSITION:");
+    G_verbose_message("XC: %.2f, YC: %.2f, ZC: %.2f", *XC, *YC, *ZC);
+    G_verbose_message("Omega %.2f, Phi %.2f, Kappa: %.2f",
+                      *Omega * 180. / M_PI, *Phi * 180. / M_PI, *Kappa * 180. / M_PI);
 
 #ifdef DEBUG
     fprintf(debug, "\nINITIAL CAMERA POSITION:\n");
@@ -282,10 +304,15 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
     epsilon.x[4][0] = *Phi;
     epsilon.x[5][0] = *Kappa;
 
+    /* adjust camera position only if enough control points are available */
+    max_iters = MAX_ITERS;
+    if (active < 4)
+	max_iters = -1;
+
 /************************** Start Iterations *****************************/
     /* itererate until convergence */
 
-    for (iter = 0; iter <= MAX_ITERS; iter++) {
+    for (iter = 0; iter < max_iters; iter++) {
 	/*  break if converged */
 	dx = delta.x[0][0];
 	dy = delta.x[1][0];
@@ -328,6 +355,7 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
 
 	/* see  WOLF 1983, Appendix, p. 591 */
 
+	/* ground -> photo */
 	M.x[0][0] = (cp * ck);
 	M.x[0][1] = (cw * sk) + (sw * sp * ck);
 	M.x[0][2] = (sw * sk) - (cw * sp * ck);
@@ -364,6 +392,12 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
 	    X = cpz->e2[i];
 	    Y = cpz->n2[i];
 	    Z = cpz->z2[i];
+	    
+	    /* adjust for earth curvature */
+	    dx = X - *XC;
+	    dy = Y - *YC;
+	    dd = (dx * dx) + (dy * dy);
+	    Z -= dd / (2.0 * ellps_a);
 
 #ifdef DEBUG
 	    fprintf(debug,
@@ -389,6 +423,31 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
 	    V = UVW.x[1][0];
 	    W = UVW.x[2][0];
 
+	    /* panorama correction
+	     * in theory either for U, V or for xbar, ybar
+	     * U, V is recommended because U, V are only used
+	     * for the residuals
+	     * correcting U, V also gives a slightly smaller RMSE */
+	    if (panorama) {
+		double a, epan;
+		
+		epan = U;
+		if (epan < 0) {
+		    a = atan2(-epan, -W);
+		    epan = -a * -W;
+		}
+		else {
+		    a = atan2(epan, -W);
+		    epan = a * -W;
+		}
+		U = epan;
+
+		V *= cos(a);
+
+		UVW.x[0][0] = U;
+		UVW.x[1][0] = V;
+	    }
+
 	    /* Form Partial derivatives of Normal Equations */
 	    xbar = x - Xp;
 	    ybar = y - Yp;
@@ -396,6 +455,7 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
 	    B.x[0][0] = (-Q1 / W) * ((xbar * m31) + (CFL * m11));
 	    B.x[0][1] = (-Q1 / W) * ((xbar * m32) + (CFL * m12));
 	    B.x[0][2] = (-Q1 / W) * ((xbar * m33) + (CFL * m13));
+
 	    B.x[0][3] = (Q1 / W) *
 		((xbar * ((-m33 * mu) + (m32 * nu))) +
 		 (CFL * ((-m13 * mu) + (m12 * nu)))
@@ -412,6 +472,7 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
 	    B.x[1][0] = (-Q1 / W) * ((ybar * m31) + (CFL * m21));
 	    B.x[1][1] = (-Q1 / W) * ((ybar * m32) + (CFL * m22));
 	    B.x[1][2] = (-Q1 / W) * ((ybar * m33) + (CFL * m23));
+
 	    B.x[1][3] = (Q1 / W) *
 		((ybar * ((-m33 * mu) + (m32 * nu))) +
 		 (CFL * ((-m23 * mu) + (m22 * nu)))
@@ -425,6 +486,7 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
 	    B.x[1][5] = (Q1 / W) *
 		(CFL * ((-m11 * lam) + (-m12 * mu) + (-m13 * nu)));
 
+	    /* residuals of image coordinates */
 	    E.x[0][0] = (-Q1) * (xbar + (CFL * (U / W)));
 	    E.x[1][0] = (-Q1) * (ybar + (CFL * (V / W)));
 
@@ -463,6 +525,8 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
 	/* Add weigth matrix of unknowns to NN */
 	m_add(&NN, &WT1, &NN);
 	/* Solve for delta */
+	
+	/* similar but not identical to solvemat in lib/imagery/georef.c */
 	if (!inverse(&NN, &N)) {	/* control point status becomes zero if matrix is non-invertable */
 	    error("Matrix Inversion (Control Points status modified)");
 	    for (i = 0; i < cpz->count; i++)
@@ -470,7 +534,6 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
 	    return (0);
 	}
 
-
 	/* delta = N-1 * CC */
 	m_mult(&N, &CC, &delta);
 	/* epsilon += delta */
@@ -484,6 +547,7 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
 #endif
 
     }  /* end ITERATION loop */
+    G_verbose_message("%d iterations to refine camera postion", iter);
 
     /* This is the solution */
     *XC = epsilon.x[0][0];
@@ -493,13 +557,16 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
     *Phi = epsilon.x[4][0];
     *Kappa = epsilon.x[5][0];
 
-    G_debug(1, "FINAL CAMERA POSITION:");
-    G_debug(1, "XC: %.2f, YC: %.2f, ZC: %.2f", *XC, *YC, *ZC);
-    G_debug(1, "Omega %.2f, Phi %.2f, Kappa: %.2f", *Omega, *Phi, *Kappa);
+    G_verbose_message("FINAL CAMERA POSITION:");
+    G_verbose_message("XC: %.2f, YC: %.2f, ZC: %.2f", *XC, *YC, *ZC);
+    G_verbose_message("Omega %.2f, Phi %.2f, Kappa: %.2f",
+                      *Omega * 180. / M_PI, *Phi * 180. / M_PI, *Kappa * 180. / M_PI);
+
     if (*ZC < 0)
 	G_warning(_("Potential BUG in ortholib: camera altitude < 0"));
 
     /*  Compute Orientation Matrix from Omega, Phi, Kappa */
+    /* ground -> photo */
     sw = sin(*Omega);
     cw = cos(*Omega);
     sp = sin(*Phi);
@@ -522,6 +589,7 @@ int I_compute_ortho_equations(struct Ortho_Control_Points *cpz,
     MO->x[2][2] = cw * cp;
 
     /* Compute Transposed Orientation Matrix from Omega, Phi, Kappa */
+    /* photo -> ground */
 
     MT->nrows = 3;
     MT->ncols = 3;
@@ -555,6 +623,7 @@ int I_ortho_ref(double e1, double n1, double z1,
     MATRIX UVW, XYZ;
     double U, V, W;
     double Xp, Yp, CFL;
+    double dx, dy, dd;
 
     /*  Initialize and zero the matrices */
     /*  Object Space Coordinates */
@@ -572,6 +641,12 @@ int I_ortho_ref(double e1, double n1, double z1,
     Yp = cam_info->Yp;
     CFL = cam_info->CFL;
 
+    /* adjust for earth curvature */
+    dx = e1 - XC;
+    dy = n1 - YC;
+    dd = (dx * dx) + (dy * dy);
+    z1 -= dd / (2.0 * ellps_a);
+
     /* Object Space (&XYZ, XC,YC,ZC, X,Y,Z); */
     XYZ.x[0][0] = e1 - XC;
     XYZ.x[1][0] = n1 - YC;
@@ -584,9 +659,29 @@ int I_ortho_ref(double e1, double n1, double z1,
     V = UVW.x[1][0];
     W = UVW.x[2][0];
 
+    /* panorama correction could also be done for e2, n2 below
+     * results are identical */
+    if (panorama) {
+	double a, epan;
+	
+	epan = U;
+	if (epan < 0) {
+	    a = atan2(-epan, -W);
+	    epan = -a * -W;
+	}
+	else {
+	    a = atan2(epan, -W);
+	    epan = a * -W;
+	}
+	U = epan;
+
+	V *= cos(a);
+    }
+
     /* This is the solution */
     *e2 = (-CFL) * (U / W);
     *n2 = (-CFL) * (V / W);
+
     return (1);
 }
 
@@ -622,6 +717,23 @@ int I_inverse_ortho_ref(double e1, double n1, double z1,
     UVW.x[0][0] = e1 - Xp;
     UVW.x[1][0] = n1 - Yp;
     UVW.x[2][0] = -CFL;
+    
+    if (panorama) {
+	double a, epan;
+	
+	epan = UVW.x[0][0];
+	if (epan < 0) {
+	    a = -epan / CFL;
+	    epan = -CFL * tan(a);
+	}
+	else {
+	    a = epan / CFL;
+	    epan = CFL * tan(a);
+	}
+	UVW.x[0][0] = epan;
+
+	UVW.x[1][0] /= cos(a);
+    }
 
     m_mult(&M, &UVW, &XYZ);