[GRASS-SVN] r44656 - in
grass/branches/releasebranch_6_4/imagery/i.ortho.photo: .
i.photo.rectify libes menu
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Dec 22 05:16:07 EST 2010
Author: mmetz
Date: 2010-12-22 02:16:07 -0800 (Wed, 22 Dec 2010)
New Revision: 44656
Added:
grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/angle.c
Modified:
grass/branches/releasebranch_6_4/imagery/i.ortho.photo/
grass/branches/releasebranch_6_4/imagery/i.ortho.photo/BUGS
grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/description.html
grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/exec.c
grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/get_wind.c
grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/local_proto.h
grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/main.c
grass/branches/releasebranch_6_4/imagery/i.ortho.photo/libes/orthoref.c
grass/branches/releasebranch_6_4/imagery/i.ortho.photo/menu/i.photo.rectify.c
Log:
merge r44551-44553,44560,44561,44569,44606,44652 from devbr6
Property changes on: grass/branches/releasebranch_6_4/imagery/i.ortho.photo
___________________________________________________________________
Modified: svn:mergeinfo
- /grass/branches/develbranch_6/imagery/i.ortho.photo:44349-44350,44472-44477,44500-44502,44507-44509,44521-44523
+ /grass/branches/develbranch_6/imagery/i.ortho.photo:44349-44350,44472-44477,44500-44502,44507-44509,44521-44523,44551-44553,44560-44561,44569,44606,44652
Modified: grass/branches/releasebranch_6_4/imagery/i.ortho.photo/BUGS
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.ortho.photo/BUGS 2010-12-22 10:10:33 UTC (rev 44655)
+++ grass/branches/releasebranch_6_4/imagery/i.ortho.photo/BUGS 2010-12-22 10:16:07 UTC (rev 44656)
@@ -13,20 +13,10 @@
window, this GCP isn't shown, should be drawn in orange color
before asking for confirmation.
-- problem in alpine regions with steep slope (or poor DEM):
-
- NOTES: the size of target window is defined in photo.rectify/global.h
- #define TIE_ROW_DIST 128
- #define TIE_COL_DIST 128
- #define NROWS 128
- #define NCOLS 128
- -> 128 is too large for alpine regions! MN 2003
- 32 works better, but slows down the processing. A dynamic selection
- based on local slope were best.
-
---------------------
TODOs
- menu 8: photo.rectify
should work without (7) photo.2target if (6) photo.init is used.
(see ERDAS)
+ requires modification of I_compute_ortho_equations() in libes/orthoref.c
Copied: grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/angle.c (from rev 44551, grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/angle.c)
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/angle.c (rev 0)
+++ grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/angle.c 2010-12-22 10:16:07 UTC (rev 44656)
@@ -0,0 +1,219 @@
+/* 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(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 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 *fbuf0, *fbuf1, *fbuf2, *tmpbuf, *outbuf;
+ int elevfd, outfd;
+ struct Cell_head cellhd;
+ struct Colors colr;
+ FCELL clr_min, clr_max;
+ struct History hist;
+ char *type;
+
+ G_message(_("Calculating camera angle to local surface..."));
+
+ select_target_env();
+
+ /* align target window to elevation map, otherwise we get artefacts
+ * like in r.slope.aspect -a */
+
+ if (G_get_cellhd(elev_name, elev_mapset, &cellhd) < 0)
+ return 0;
+
+ G_align_window(&target_window, &cellhd);
+ G_set_window(&target_window);
+
+ elevfd = G_open_cell_old(elev_name, elev_mapset);
+ if (elevfd < 0) {
+ G_fatal_error(_("Could not open elevation raster"));
+ return 1;
+ }
+
+ nrows = target_window.rows;
+ ncols = target_window.cols;
+
+ outfd = G_open_raster_new(name, FCELL_TYPE);
+ fbuf0 = G_allocate_raster_buf(FCELL_TYPE);
+ fbuf1 = G_allocate_raster_buf(FCELL_TYPE);
+ fbuf2 = G_allocate_raster_buf(FCELL_TYPE);
+ outbuf = 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;
+ G_get_raster_row(elevfd, fbuf1, 0, FCELL_TYPE);
+ G_get_raster_row(elevfd, fbuf2, 1, FCELL_TYPE);
+
+ for (row = 0; row < nrows; row++) {
+ G_percent(row, nrows, 2);
+
+ G_set_null_value(outbuf, ncols, FCELL_TYPE);
+
+ /* first and last row */
+ if (row == 0 || row == nrows - 1) {
+ G_put_raster_row(outfd, outbuf, FCELL_TYPE);
+ continue;
+ }
+
+ tmpbuf = fbuf0;
+ fbuf0 = fbuf1;
+ fbuf1 = fbuf2;
+ fbuf2 = tmpbuf;
+
+ G_get_raster_row(elevfd, fbuf2, row + 1, FCELL_TYPE);
+
+ north = G_row_to_northing(row + 0.5, &target_window);
+
+ for (col = 1; col < ncols - 1; col++) {
+
+ e1 = fbuf0[col - 1];
+ if (G_is_d_null_value(&e1))
+ continue;
+ e2 = fbuf0[col];
+ if (G_is_d_null_value(&e2))
+ continue;
+ e3 = fbuf0[col + 1];
+ if (G_is_d_null_value(&e3))
+ continue;
+ e4 = fbuf1[col - 1];
+ if (G_is_d_null_value(&e4))
+ continue;
+ e5 = fbuf1[col];
+ if (G_is_d_null_value(&e5))
+ continue;
+ e6 = fbuf1[col + 1];
+ if (G_is_d_null_value(&e6))
+ continue;
+ e7 = fbuf2[col - 1];
+ if (G_is_d_null_value(&e7))
+ continue;
+ e8 = fbuf2[col];
+ if (G_is_d_null_value(&e8))
+ continue;
+ e9 = fbuf2[col + 1];
+ if (G_is_d_null_value(&e9))
+ continue;
+
+ 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));
+
+ outbuf[col] = c_angle * radians_to_degrees;
+ if (c_angle_min > outbuf[col])
+ c_angle_min = outbuf[col];
+ }
+ G_put_raster_row(outfd, outbuf, FCELL_TYPE);
+ }
+ G_percent(row, nrows, 2);
+
+ G_close_cell(elevfd);
+ G_close_cell(outfd);
+ G_free(fbuf0);
+ G_free(fbuf1);
+ G_free(fbuf2);
+ G_free(outbuf);
+
+ type = "raster";
+ G_short_history(name, type, &hist);
+ G_command_history(&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;
+}
Modified: grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/description.html
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/description.html 2010-12-22 10:10:33 UTC (rev 44655)
+++ grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/description.html 2010-12-22 10:16:07 UTC (rev 44656)
@@ -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/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/exec.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/exec.c 2010-12-22 10:10:33 UTC (rev 44655)
+++ grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/exec.c 2010-12-22 10:16:07 UTC (rev 44656)
@@ -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;
@@ -100,6 +100,7 @@
/* Initialze History */
type = "raster";
G_short_history(name, type, &hist);
+ G_command_history(&hist);
G_write_history(result, &hist);
select_current_env();
@@ -111,9 +112,14 @@
G_free(result);
}
+
close(ebuffer->fd);
release_cache(ebuffer);
+ if (angle_map) {
+ camera_angle(angle_map);
+ }
+
G_done_msg(" ");
return 0;
Modified: grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/get_wind.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/get_wind.c 2010-12-22 10:10:33 UTC (rev 44655)
+++ grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/get_wind.c 2010-12-22 10:16:07 UTC (rev 44656)
@@ -62,6 +62,9 @@
} nw, ne, se, sw;
/* get an average elevation from the active control points */
+ /* for mountainous areas this is a very rough approximation
+ * which would become more accurate only if actual elevation
+ * values are used */
get_aver_elev(&group.control_points, &aver_z);
G_debug(1, "Aver elev = %f", aver_z);
@@ -94,7 +97,6 @@
w1->north, e0, n0);
G_debug(1, "target x = %f y = %f", e, n);
-
ne.n = n;
ne.e = e;
if (n > w2->north)
Modified: grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/local_proto.h
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/local_proto.h 2010-12-22 10:10:33 UTC (rev 44655)
+++ grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/local_proto.h 2010-12-22 10:16:07 UTC (rev 44656)
@@ -1,3 +1,6 @@
+/* angle.c */
+int camera_angle(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/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/main.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/main.c 2010-12-22 10:10:33 UTC (rev 44655)
+++ grass/branches/releasebranch_6_4/imagery/i.ortho.photo/i.photo.rectify/main.c 2010-12-22 10:16:07 UTC (rev 44656)
@@ -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;
@@ -98,7 +100,7 @@
ext->description = _("Output raster map(s) suffix");
tres = G_define_option();
- tres->key = "res";
+ tres->key = "resolution";
tres->type = TYPE_DOUBLE;
tres->required = NO;
tres->description = _("Target resolution (ignored if -c flag used)");
@@ -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);
}
Modified: grass/branches/releasebranch_6_4/imagery/i.ortho.photo/libes/orthoref.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.ortho.photo/libes/orthoref.c 2010-12-22 10:10:33 UTC (rev 44655)
+++ grass/branches/releasebranch_6_4/imagery/i.ortho.photo/libes/orthoref.c 2010-12-22 10:16:07 UTC (rev 44656)
@@ -30,8 +30,9 @@
#include <signal.h>
#include <stdio.h>
#include <math.h>
-#include "orthophoto.h"
#include <grass/imagery.h>
+#include <grass/glocale.h>
+#include "orthophoto.h"
#include "matrixdefs.h"
#include "local_proto.h"
@@ -54,9 +55,6 @@
MATRIX *MO, MATRIX *MT)
{
MATRIX delta, epsilon, B, BT, C, E, N, CC, NN, UVW, XYZ, M, WT1;
- double meanx, meany;
- double X1 = 0, X2 = 0, x1 = 0, x2 = 0, Z1 = 0, Y1 = 0, Y2 = 0, y1 =
- 0, y2 = 0, dist_grnd = 0, dist_photo = 0;
double x, y, z, X, Y, Z, Xp, Yp, CFL;
double lam, mu, nu, U, V, W;
double xbar, ybar;
@@ -187,14 +185,13 @@
WT1.x[5][5] = (Q1 / (init_info->kappa_var * init_info->kappa_var));
}
else { /* set from mean values of active control points */
- double meanX, meanY, meanZ;
- int nz;
-
+ double dist_grnd, dist_photo;
+ double meanx, meany, meanX, meanY, meanZ;
+
/* set intial XC and YC from mean values of control points */
- meanx = meany = 0;
- meanX = meanY = meanZ = 0;
- *ZC = 0;
- n = nz = 0;
+ meanx = meany = meanX = meanY = meanZ = 0;
+ x = y = X = Y = 0;
+ n = 0;
first = 1;
for (i = 0; i < cpz->count; i++) {
if (cpz->status[i] <= 0)
@@ -202,71 +199,52 @@
/* set initial XC, YC */
if (first) {
- X1 = *(cpz->e2);
- x1 = *(cpz->e1);
- Y1 = *(cpz->n2);
- y1 = *(cpz->n1);
- Z1 = *(cpz->z2);
+ X = cpz->e2[i];
+ x = cpz->e1[i];
+ Y = cpz->n2[i];
+ y = cpz->n1[i];
first = 0;
}
- if (!first) {
- X2 = *(cpz->e2);
- x2 = *(cpz->e1);
- Y2 = *(cpz->n2);
- y2 = *(cpz->n1);
- Z1 = *(cpz->z2);
- dist_photo =
- sqrt(((x1 - x2) * (x1 - x2)) + ((y1 - y2) * (y1 - y2)));
- dist_grnd =
- sqrt(((X1 - X2) * (X1 - X2)) + ((Y1 - Y2) * (Y1 - Y2)));
+ else {
+ /* set initial ZC from: */
+ /* scale ~= dist_photo/dist_grnd ~= (CFL)/(Z - ZC) */
+ /* ZC ~= Z + CFL(dist_grnd)/(dist_photo) */
+ dx = cpz->e1[i] - x;
+ dy = cpz->n1[i] - y;
+ dist_photo = sqrt(dx * dx + dy * dy);
+ dx = cpz->e2[i] - X;
+ dy = cpz->n2[i] - Y;
+ dist_grnd = sqrt(dx * dx + dy * dy);
+
+ if (dist_photo != 0 && dist_grnd != 0) {
+ meanZ += cpz->z2[i] + (CFL * (dist_grnd) / (dist_photo));
+ meanx += cpz->e1[i];
+ meany += cpz->n1[i];
+ meanX += cpz->e2[i];
+ meanY += cpz->n2[i];
- /* set initial ZC from: */
- /* scale ~= dist_photo/dist_grnd ~= (CFL)/(Z1 - Zc) */
- /* Zc ~= Z1 + CFL(dist_grnd)/(dist_photo) */
-
- if (dist_grnd != 0 && dist_photo != 0) {
- meanZ += Z1 + (CFL * (dist_grnd) / (dist_photo));
- nz++;
+ n++;
}
}
-
- n++;
- meanx += *((cpz->e1)++);
- meany += *((cpz->n1)++);
- meanX += *((cpz->e2)++);
- meanY += *((cpz->n2)++);
- (cpz->z2)++;
}
- *XC = meanX / n;
- *YC = meanY / n;
- *ZC = meanZ / (nz - 1);
- G_debug(2, "XC: %.2f, YC: %.2f, ZC: %.2f", *XC, *YC, *ZC);
- meanx = (meanx - x1) / (n - 1);
- meany = (meany - y1) / (n - 1);
- meanX = (meanX - X1) / (n - 1);
- meanY = (meany - Y1) / (n - 1);
+ if (!n) /* Poorly placed Control Points */
+ return -1;
- /* reset pointers */
- for (i = 0; i < cpz->count; i++) {
- if (cpz->status[i] <= 0)
- continue;
- (cpz->e1)--;
- (cpz->e2)--;
- (cpz->n1)--;
- (cpz->n2)--;
- (cpz->z2)--;
- }
+ meanx /= n;
+ meany /= n;
+ meanX /= n;
+ meanY /= n;
+ meanZ /= n;
+ *XC = meanX;
+ *YC = meanY;
+ *ZC = meanZ;
/* set initial rotations to zero (radians) */
*Omega = *Phi = 0.0;
/* find an initial kappa */
- /*
- kappa1 = atan2((y2 - y1), (x2 - x1));
- kappa2 = atan2((Y2 - Y1), (X2 - X1));
- */
- kappa1 = atan2((meany - y1), (meanx - x1));
- kappa2 = atan2((meanY - Y1), (meanX - X1));
+ kappa1 = atan2((meany - y), (meanx - x));
+ kappa2 = atan2((meanY - Y), (meanX - X));
*Kappa = (kappa2 - kappa1);
/* set initial weights */
@@ -284,6 +262,9 @@
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);
#ifdef DEBUG
fprintf(debug, "\nINITIAL CAMERA POSITION:\n");
@@ -373,29 +354,28 @@
fprintf(debug, "\n\t\t\tIn Summation count = %d \n", i);
#endif
- x = *((cpz->e1)++);
- y = *((cpz->n1)++);
- z = *((cpz->z1)++);
- X = *((cpz->e2)++);
- Y = *((cpz->n2)++);
- Z = *((cpz->z2)++);
-
if (cpz->status[i] <= 0)
continue;
+ x = cpz->e1[i];
+ y = cpz->n1[i];
+ z = cpz->z1[i];
+ X = cpz->e2[i];
+ Y = cpz->n2[i];
+ Z = cpz->z2[i];
+
#ifdef DEBUG
fprintf(debug,
"\n\t\t\timage:\n \t\t\tx = \t%f \n\t\t\ty = \t%f \n\t\t\tz = \t%f \n\t\t\tobject:\n \t\t\tX = \t%f \n\t\t\tY = \t%f \n\t\t\tZ = \t%f \n",
x, y, z, X, Y, Z);
#endif
-
/* Compute Obj. Space coordinates */
XYZ.x[0][0] = X - *XC;
XYZ.x[1][0] = Y - *YC;
XYZ.x[2][0] = Z - *ZC;
- /* just an abbreviations */
+ /* just an abbreviation */
lam = XYZ.x[0][0];
mu = XYZ.x[1][0];
nu = XYZ.x[2][0];
@@ -466,16 +446,6 @@
m_add(&CC, &C, &CC);
} /* end summation loop over all active control points */
- /* reset pointers */
- for (i = 0; i < cpz->count; i++) {
- (cpz->e1)--;
- (cpz->n1)--;
- (cpz->z1)--;
- (cpz->e2)--;
- (cpz->n2)--;
- (cpz->z2)--;
- }
-
#ifdef DEBUG
fprintf(debug, "\n\tN: \n");
for (i = 0; i < 6; i++)
@@ -522,7 +492,13 @@
*Phi = epsilon.x[4][0];
*Kappa = epsilon.x[5][0];
- /* Compute Orientation Matrix M from (Omega, Phi, Kappa); */
+ 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);
+ if (*ZC < 0)
+ G_warning(_("Potential BUG in ortholib: camera altitude < 0"));
+
+ /* Compute Orientation Matrix from Omega, Phi, Kappa */
sw = sin(*Omega);
cw = cos(*Omega);
sp = sin(*Phi);
@@ -544,13 +520,13 @@
MO->x[2][1] = -(sw * cp);
MO->x[2][2] = cw * cp;
- /* Compute Transposed Orientation Matrix (Omega, Phi, Kappa); */
+ /* Compute Transposed Orientation Matrix from Omega, Phi, Kappa */
MT->nrows = 3;
MT->ncols = 3;
zero(MT);
- /* M Transposed */
+ /* Transposed Matrix */
MT->x[0][0] = cp * ck;
MT->x[1][0] = cw * sk + (sw * sp * ck);
MT->x[2][0] = sw * sk - (cw * sp * ck);
@@ -590,12 +566,12 @@
zero(&UVW);
/************************ Start the work ******************************/
- /* set Xp, Yp, and CFL from cam_info */
+ /* Set Xp, Yp, and CFL from cam_info */
Xp = cam_info->Xp;
Yp = cam_info->Yp;
CFL = cam_info->CFL;
- /* ObjSpace (&XYZ, XC,YC,ZC, X,Y,Z); */
+ /* Object Space (&XYZ, XC,YC,ZC, X,Y,Z); */
XYZ.x[0][0] = e1 - XC;
XYZ.x[1][0] = n1 - YC;
XYZ.x[2][0] = z1 - ZC;
@@ -636,19 +612,19 @@
zero(&UVW);
/********************** Start the work ********************************/
- /* set Xp, Yp, and CFL from cam_info */
+ /* Set Xp, Yp, and CFL from cam_info */
Xp = cam_info->Xp;
Yp = cam_info->Yp;
CFL = cam_info->CFL;
- /* ImageSpace */
+ /* Image Space */
UVW.x[0][0] = e1 - Xp;
UVW.x[1][0] = n1 - Yp;
UVW.x[2][0] = -CFL;
m_mult(&M, &UVW, &XYZ);
- /* Image Space */
+ /* Object Space */
lam = XYZ.x[0][0];
mu = XYZ.x[1][0];
nu = XYZ.x[2][0];
Modified: grass/branches/releasebranch_6_4/imagery/i.ortho.photo/menu/i.photo.rectify.c
===================================================================
--- grass/branches/releasebranch_6_4/imagery/i.ortho.photo/menu/i.photo.rectify.c 2010-12-22 10:10:33 UTC (rev 44655)
+++ grass/branches/releasebranch_6_4/imagery/i.ortho.photo/menu/i.photo.rectify.c 2010-12-22 10:16:07 UTC (rev 44656)
@@ -117,6 +117,40 @@
}
sprintf(pgm, "%s extension=%s", pgm, extension);
+ /* compute local camera angle */
+ G_clear_screen();
+ ok = G_yes(_("\nCompute local camera angle? "), 0);
+ if (ok) {
+ char angle_name[GNAME_MAX];
+
+ sprintf(angle_name, "%s.camera_angle", groupname);
+
+ repeat = 1;
+ while (repeat) {
+ repeat = 0;
+ V_clear();
+ V_line(1, _("Enter a name for the camera angle map:"));
+ V_ques(angle_name, 's', 3, 0, 30);
+ V_intrpt_ok();
+ if (!V_call())
+ return 0;
+
+ /* test for legal file name */
+ if (G_legal_filename(angle_name) < 0) {
+ G_clear_screen();
+ fprintf(stderr, _("Map name <%s> is illegal"), angle_name);
+ repeat = G_yes(_("\nChoose another name? "), 1);
+ if (!repeat) {
+ fprintf(stderr,_("Orthorectification cancelled."));
+ G_sleep(3);
+ return 0;
+ }
+ }
+ else
+ sprintf(pgm, "%s angle=%s", pgm, angle_name);
+ }
+ }
+
/* overwrite maps in target location/mapset */
G_clear_screen();
ok = G_yes(_("\nOverwrite maps in target location/mapset? "), 0);
@@ -157,9 +191,9 @@
while (1) {
char buf[100];
- fprintf(stderr, "Desired target resolution\n");
+ fprintf(stderr, "Enter desired target resolution, or\n");
fprintf(stderr,
- " RETURN determine automatically\n");
+ " RETURN to determine it automatically:\n");
fprintf(stderr, "> ");
if (!G_gets(buf))
continue;
@@ -254,7 +288,7 @@
char buf[100];
int seg_mb = max_mb + 0.5;
- fprintf(stderr, _("Amount of memory to use in MB\n"));
+ fprintf(stderr, _("Enter amount of memory to use in MB, or\n"));
if (use_target_window)
fprintf(stderr, _("RETURN use %d MB to keep all data in RAM\n"), seg_mb);
else {
More information about the grass-commit
mailing list