[GRASS-SVN] r44551 - grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify

svn_grass at osgeo.org svn_grass at osgeo.org
Tue Dec 7 03:23:03 EST 2010


Author: mmetz
Date: 2010-12-07 00:23:03 -0800 (Tue, 07 Dec 2010)
New Revision: 44551

Added:
   grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/angle.c
Modified:
   grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/description.html
   grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/exec.c
   grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/local_proto.h
   grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/main.c
Log:
add new option for camera angle

Added: grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/angle.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/angle.c	                        (rev 0)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/angle.c	2010-12-07 08:23:03 UTC (rev 44551)
@@ -0,0 +1,195 @@
+/* calculate camera angle to local surface
+ * 90 degrees: orthogonal to local surface
+ * 0 degress: parallel to local surface
+ * < 0 degrees: not visible by camera
+ * 
+ * Earth curvature is not considered, assuming that the extends of the 
+ * imagery to be orthorectified are rather small
+ * Shadowing effects by ridges and peaks are not considered */
+
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+#include <math.h>
+#include "global.h"
+
+int camera_angle(struct cache *ebuffer, char *name)
+{
+    int row, col, nrows, ncols;
+    double XC = group.XC;
+    double YC = group.YC;
+    double ZC = group.ZC;
+    double c_angle, c_angle_min, c_alt, c_az, slope, aspect;
+    double radians_to_degrees = 180.0 / M_PI;
+    /* double degrees_to_radians = M_PI / 180.0; */
+    DCELL *ep;
+    DCELL e1, e2, e3, e4, e5, e6, e7, e8, e9;
+    double factor, V, H, dx, dy, dz, key;
+    double north, south, east, west, ns_med;
+    FCELL *fbuf;
+    int outfd;
+    struct Colors colr;
+    FCELL clr_min, clr_max;
+    struct History hist;
+    char *type;
+
+    G_message(_("Calculating camera angle to local surface..."));
+    
+    select_target_env();
+    
+    nrows = target_window.rows;
+    ncols = target_window.cols;
+    
+    outfd = G_open_raster_new(name, FCELL_TYPE);
+    fbuf = G_allocate_raster_buf(FCELL_TYPE);
+    
+    /* give warning if location units are different from meters and zfactor=1 */
+    factor = G_database_units_to_meters_factor();
+    if (factor != 1.0)
+	G_warning(_("Converting units to meters, factor=%.6f"), factor);
+
+    G_begin_distance_calculations();
+    north = G_row_to_northing(0.5, &target_window);
+    ns_med = G_row_to_northing(1.5, &target_window);
+    south = G_row_to_northing(2.5, &target_window);
+    east = G_col_to_easting(2.5, &target_window);
+    west = G_col_to_easting(0.5, &target_window);
+    V = G_distance(east, north, east, south) * 4;
+    H = G_distance(east, ns_med, west, ns_med) * 4;
+    
+    c_angle_min = 90;
+    for (row = 0; row < nrows; row++) {
+	G_percent(row, nrows, 2);
+	
+	G_set_null_value(fbuf, ncols, FCELL_TYPE);
+
+	if (row == 0 || row == nrows - 1) {
+	    G_put_raster_row(outfd, fbuf, FCELL_TYPE);
+	    continue;
+	}
+
+	north = G_row_to_northing(row + 0.5, &target_window);
+
+	for (col = 1; col < ncols - 1; col++) {
+	    
+	    ep = CPTR(ebuffer, row - 1, col - 1);
+	    if (G_is_d_null_value(ep))
+		continue;
+	    e1 = *ep;
+	    ep = CPTR(ebuffer, row - 1, col);
+	    if (G_is_d_null_value(ep))
+		continue;
+	    e2 = *ep;
+	    ep = CPTR(ebuffer, row - 1, col + 1);
+	    if (G_is_d_null_value(ep))
+		continue;
+	    e3 = *ep;
+	    ep = CPTR(ebuffer, row, col - 1);
+	    if (G_is_d_null_value(ep))
+		continue;
+	    e4 = *ep;
+	    ep = CPTR(ebuffer, row, col);
+	    if (G_is_d_null_value(ep))
+		continue;
+	    e5 = *ep;
+	    ep = CPTR(ebuffer, row, col + 1);
+	    if (G_is_d_null_value(ep))
+		continue;
+	    e6 = *ep;
+	    ep = CPTR(ebuffer, row + 1, col - 1);
+	    if (G_is_d_null_value(ep))
+		continue;
+	    e7 = *ep;
+	    ep = CPTR(ebuffer, row + 1, col);
+	    if (G_is_d_null_value(ep))
+		continue;
+	    e8 = *ep;
+	    ep = CPTR(ebuffer, row + 1, col + 1);
+	    if (G_is_d_null_value(ep))
+		continue;
+	    e9 = *ep;
+	    
+	    dx = ((e1 + e4 + e4 + e7) - (e3 + e6 + e6 + e9)) / H;
+	    dy = ((e7 + e8 + e8 + e9) - (e1 + e2 + e2 + e3)) / V;
+	    
+	    /* compute topographic parameters */
+	    key = dx * dx + dy * dy;
+	    /* slope in radians */
+	    slope = atan(sqrt(key));
+
+	    /* aspect in radians */
+	    if (key == 0.)
+		aspect = 0.;
+	    else if (dx == 0) {
+		if (dy > 0)
+		    aspect = M_PI / 2;
+		else
+		    aspect = 1.5 * M_PI;
+	    }
+	    else {
+		aspect = atan2(dy, dx);
+		if (aspect <= 0.)
+		    aspect = 2 * M_PI + aspect;
+	    }
+	    
+	    /* camera altitude angle in radians */
+	    east = G_col_to_easting(col + 0.5, &target_window);
+	    dx = east - XC;
+	    dy = north - YC;
+	    dz = ZC - e5;
+	    c_alt = atan(sqrt(dx * dx + dy * dy) / dz);
+
+	    /* camera azimuth angle in radians */
+	    c_az = atan(dy / dx);
+	    if (east < XC && north != YC)
+		c_az += M_PI;
+	    else if (north < YC && east > XC)
+		c_az += 2 * M_PI;
+		
+	    /* camera angle to real ground */
+	    /* orthogonal to ground: 90 degrees */
+	    /* parallel to ground: 0 degrees */
+	    c_angle = asin(cos(c_alt) * cos(slope) - sin(c_alt) * sin(slope) * cos(c_az - aspect));
+	    
+	    fbuf[col] = c_angle * radians_to_degrees;
+	    if (c_angle_min > fbuf[col])
+		c_angle_min = fbuf[col];
+	}
+	G_put_raster_row(outfd, fbuf, FCELL_TYPE);
+    }
+    G_percent(row, nrows, 2);
+
+    G_close_cell(outfd);
+    G_free(fbuf);
+
+    type = "raster";
+    G_short_history(name, type, &hist);
+    G_write_history(name, &hist);
+    
+    G_init_colors(&colr);
+    if (c_angle_min < 0) {
+	clr_min = (FCELL)((int)(c_angle_min / 10 - 1)) * 10;
+	clr_max = 0;
+	G_add_f_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 0,
+				  0, 0, &colr);
+    }
+    clr_min = 0;
+    clr_max = 10;
+    G_add_f_raster_color_rule(&clr_min, 0, 0, 0, &clr_max, 255,
+			      0, 0, &colr);
+    clr_min = 10;
+    clr_max = 40;
+    G_add_f_raster_color_rule(&clr_min, 255, 0, 0, &clr_max, 255,
+			      255, 0, &colr);
+    clr_min = 40;
+    clr_max = 90;
+    G_add_f_raster_color_rule(&clr_min, 255, 255, 0, &clr_max, 0,
+			      255, 0, &colr);
+
+    G_write_colors(name, G_mapset(), &colr);
+
+
+    select_current_env();
+
+    return 1;
+}


Property changes on: grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/angle.c
___________________________________________________________________
Added: svn:mime-type
   + text/x-csrc
Added: svn:eol-style
   + native

Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/description.html
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/description.html	2010-12-06 21:11:31 UTC (rev 44550)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/description.html	2010-12-07 08:23:03 UTC (rev 44551)
@@ -27,6 +27,17 @@
 The rectified image will be located in the target LOCATION when the program
 is completed. The original unrectified files are not modified or removed.
 <P>
+The optional <em>angle</em> output holds the camera angle in degrees to
+the local surface, considering local slope and aspect. A value of 90 
+degrees indicates that the camera angle was orthogonal to the local 
+surface, a value of 0 degrees indicates that the camera angle was 
+parallel to the local surface and negative values indicate that the 
+surface was invisible to the camera. As a rule of thumb, values below 30
+degrees indicate problem areas where the orthorectified output will 
+appear blurred. Because terrain shadowing effects are not considered, 
+areas with high camera angles may also appear blurred if they are located
+(viewed from the camera position) behind mountain ridges or peaks.
+<P>
 <EM>i.photo.rectify</EM> can be run directly, specifying options in the 
 command line or the GUI, or it can be invoked as OPTION 8 through 
 <A HREF="i.ortho.photo.html">i.ortho.photo</A>. If invoked though 
@@ -38,10 +49,13 @@
 be rectified. If this option is not chosen, you are asked to specify for 
 each image within the imagery group whether it should be rectified or not.
 <P>
-More than one file may be rectified at a time.  Each file
+More than one file may be rectified at a time. Each file
 should have a unique output file name. The next prompt asks you for an 
 extension to be appended to the rectified images.
 <P>
+The next prompt will ask you whether a camera angle map should be 
+produced and if yes, what should be its name.
+<P>
 After that you are asked if overwriting existing maps in the target 
 location and mapset should be allowed.
 <P>

Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/exec.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/exec.c	2010-12-06 21:11:31 UTC (rev 44550)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/exec.c	2010-12-07 08:23:03 UTC (rev 44551)
@@ -13,7 +13,7 @@
 #include <fcntl.h>
 #include "global.h"
 
-int exec_rectify(char *extension, char *interp_method)
+int exec_rectify(char *extension, char *interp_method, char *angle_map)
 {
     char *name;
     char *mapset;
@@ -111,6 +111,11 @@
 
 	G_free(result);
     }
+    
+    if (angle_map) {
+	camera_angle(ebuffer, angle_map);
+    }
+    
     close(ebuffer->fd);
     release_cache(ebuffer);
 

Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/local_proto.h
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/local_proto.h	2010-12-06 21:11:31 UTC (rev 44550)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/local_proto.h	2010-12-07 08:23:03 UTC (rev 44551)
@@ -1,3 +1,6 @@
+/* angle.c */
+int camera_angle(struct cache *, char *);
+
 /* aver_z.c */
 int get_aver_elev(struct Ortho_Control_Points *, double *);
 
@@ -18,7 +21,7 @@
 int Compute_ref_equation(void);
 
 /* exec.c */
-int exec_rectify(char *, char *);
+int exec_rectify(char *, char *, char *);
 
 /* get_wind.c */
 int get_ref_window(struct Cell_head *);

Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/main.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/main.c	2010-12-06 21:11:31 UTC (rev 44550)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/main.c	2010-12-07 08:23:03 UTC (rev 44551)
@@ -73,8 +73,10 @@
      *ext,			/* extension */
      *tres,			/* target resolution */
      *mem,			/* amount of memory for cache */
-     *interpol;			/* interpolation method:
+     *interpol,			/* interpolation method:
 				   nearest neighbor, bilinear, cubic */
+     *angle;			/* camera angle relative to ground surface */
+
     struct Flag *c, *a;
     struct GModule *module;
 
@@ -120,6 +122,11 @@
     interpol->answer = "nearest";
     interpol->options = ipolname;
     interpol->description = _("Interpolation method to use");
+    
+    angle = G_define_standard_option(G_OPT_R_OUTPUT);
+    angle->key = "angle";
+    angle->required = NO;
+    angle->description = _("Raster map with camera angle relative to ground surface");
 
     c = G_define_flag();
     c->key = 'c';
@@ -187,25 +194,39 @@
 	}
     }
     else {
-	char xname[GNAME_MAX], xmapset[GMAPSET_MAX], *name;
+	char xname[GNAME_MAX], xmapset[GMAPSET_MAX], *name, *mapset;
 
 	for (n = 0; n < group.group_ref.nfiles; n++)
 		ref_list[n] = 0;
 
 	for (m = 0; m < k; m++) {
 	    got_file = 0;
-	    if (G__name_is_fully_qualified(ifile->answers[m], xname, xmapset))
+	    if (G__name_is_fully_qualified(ifile->answers[m], xname, xmapset)) {
 		name = xname;
-	    else
+		mapset = xmapset;
+	    }
+	    else {
 		name = ifile->answers[m];
+		mapset = NULL;
+	    }
 
 	    got_file = 0;
 	    for (n = 0; n < group.group_ref.nfiles; n++) {
-		if (strcmp(name, group.group_ref.file[n].name) == 0) {
-		    got_file = 1;
-		    ref_list[n] = 1;
-		    break;
+		if (mapset) {
+		    if (strcmp(name, group.group_ref.file[n].name) == 0 &&
+		        strcmp(mapset, group.group_ref.file[n].mapset) == 0) {
+			got_file = 1;
+			ref_list[n] = 1;
+			break;
+		    }
 		}
+		else {
+		    if (strcmp(name, group.group_ref.file[n].name) == 0) {
+			got_file = 1;
+			ref_list[n] = 1;
+			break;
+		    }
+		}
 	    }
 	    if (got_file == 0)
 		err_exit(ifile->answers[m], group.name);
@@ -264,6 +285,16 @@
 		G_fatal_error(_("Orthorectification cancelled."));
 	    }
 	}
+	if (angle->answer) {
+	    if (G_find_cell(angle->answer, G_mapset())) {
+		G_warning(_("The following raster map already exists in"));
+		G_warning(_("target LOCATION %s, MAPSET %s:"),
+			  G_location(), G_mapset());
+		G_warning("<%s>", angle->answer);
+		G_fatal_error(_("Orthorectification cancelled."));
+	    }
+	}
+	
 	select_current_env();
     }
     else
@@ -337,7 +368,7 @@
     }
 
     /* go do it */
-    exec_rectify(extension, interpol->answer);
+    exec_rectify(extension, interpol->answer, angle->answer);
 
     exit(EXIT_SUCCESS);
 }



More information about the grass-commit mailing list