[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