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

@@ -2242,7 +2242,7 @@ def _Create(self):
                         _('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 @@ def _Create(self):
                         _('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 @@ <h2>DESCRIPTION</h2>
 	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 @@ <h2>DESCRIPTION</h2>
 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 @@ <h2>EXAMPLE</h2>
     <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 @@ <h2>EXAMPLE</h2>
 </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;
 
 
@@ -134,6 +134,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);
 
@@ -251,15 +255,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);