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

svn_grass at osgeo.org svn_grass at osgeo.org
Tue Dec 7 12:49:57 EST 2010


Author: mmetz
Date: 2010-12-07 09:49:57 -0800 (Tue, 07 Dec 2010)
New Revision: 44560

Modified:
   grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/angle.c
   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
Log:
always align camera angle map to elevation map to avoid artifacts

Modified: 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	2010-12-07 14:19:08 UTC (rev 44559)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/angle.c	2010-12-07 17:49:57 UTC (rev 44560)
@@ -13,7 +13,7 @@
 #include <math.h>
 #include "global.h"
 
-int camera_angle(struct cache *ebuffer, char *name)
+int camera_angle(char *name)
 {
     int row, col, nrows, ncols;
     double XC = group.XC;
@@ -22,12 +22,12 @@
     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;
+    FCELL *fbuf0, *fbuf1, *fbuf2, *tmpbuf, *outbuf;
+    int elevfd, outfd;
+    struct Cell_head cellhd;
     struct Colors colr;
     FCELL clr_min, clr_max;
     struct History hist;
@@ -37,11 +37,29 @@
     
     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);
-    fbuf = G_allocate_raster_buf(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();
@@ -58,56 +76,58 @@
     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(fbuf, ncols, FCELL_TYPE);
+	G_set_null_value(outbuf, ncols, FCELL_TYPE);
 
+	/* first and last row */
 	if (row == 0 || row == nrows - 1) {
-	    G_put_raster_row(outfd, fbuf, FCELL_TYPE);
+	    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++) {
 	    
-	    ep = CPTR(ebuffer, row - 1, col - 1);
-	    if (G_is_d_null_value(ep))
+	    e1 = fbuf0[col - 1];
+	    if (G_is_d_null_value(&e1))
 		continue;
-	    e1 = *ep;
-	    ep = CPTR(ebuffer, row - 1, col);
-	    if (G_is_d_null_value(ep))
+	    e2 = fbuf0[col];
+	    if (G_is_d_null_value(&e2))
 		continue;
-	    e2 = *ep;
-	    ep = CPTR(ebuffer, row - 1, col + 1);
-	    if (G_is_d_null_value(ep))
+	    e3 = fbuf0[col + 1];
+	    if (G_is_d_null_value(&e3))
 		continue;
-	    e3 = *ep;
-	    ep = CPTR(ebuffer, row, col - 1);
-	    if (G_is_d_null_value(ep))
+	    e4 = fbuf1[col - 1];
+	    if (G_is_d_null_value(&e4))
 		continue;
-	    e4 = *ep;
-	    ep = CPTR(ebuffer, row, col);
-	    if (G_is_d_null_value(ep))
+	    e5 = fbuf1[col];
+	    if (G_is_d_null_value(&e5))
 		continue;
-	    e5 = *ep;
-	    ep = CPTR(ebuffer, row, col + 1);
-	    if (G_is_d_null_value(ep))
+	    e6 = fbuf1[col + 1];
+	    if (G_is_d_null_value(&e6))
 		continue;
-	    e6 = *ep;
-	    ep = CPTR(ebuffer, row + 1, col - 1);
-	    if (G_is_d_null_value(ep))
+	    e7 = fbuf2[col - 1];
+	    if (G_is_d_null_value(&e7))
 		continue;
-	    e7 = *ep;
-	    ep = CPTR(ebuffer, row + 1, col);
-	    if (G_is_d_null_value(ep))
+	    e8 = fbuf2[col];
+	    if (G_is_d_null_value(&e8))
 		continue;
-	    e8 = *ep;
-	    ep = CPTR(ebuffer, row + 1, col + 1);
-	    if (G_is_d_null_value(ep))
+	    e9 = fbuf2[col + 1];
+	    if (G_is_d_null_value(&e9))
 		continue;
-	    e9 = *ep;
 	    
 	    dx = ((e1 + e4 + e4 + e7) - (e3 + e6 + e6 + e9)) / H;
 	    dy = ((e7 + e8 + e8 + e9) - (e1 + e2 + e2 + e3)) / V;
@@ -151,16 +171,20 @@
 	    /* 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];
+	    outbuf[col] = c_angle * radians_to_degrees;
+	    if (c_angle_min > outbuf[col])
+		c_angle_min = outbuf[col];
 	}
-	G_put_raster_row(outfd, fbuf, FCELL_TYPE);
+	G_put_raster_row(outfd, outbuf, FCELL_TYPE);
     }
     G_percent(row, nrows, 2);
 
+    G_close_cell(elevfd);
     G_close_cell(outfd);
-    G_free(fbuf);
+    G_free(fbuf0);
+    G_free(fbuf1);
+    G_free(fbuf2);
+    G_free(outbuf);
 
     type = "raster";
     G_short_history(name, type, &hist);
@@ -188,7 +212,6 @@
 
     G_write_colors(name, G_mapset(), &colr);
 
-
     select_current_env();
 
     return 1;

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-07 14:19:08 UTC (rev 44559)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/exec.c	2010-12-07 17:49:57 UTC (rev 44560)
@@ -112,13 +112,13 @@
 	G_free(result);
     }
     
+    close(ebuffer->fd);
+    release_cache(ebuffer);
+
     if (angle_map) {
-	camera_angle(ebuffer, angle_map);
+	camera_angle(angle_map);
     }
     
-    close(ebuffer->fd);
-    release_cache(ebuffer);
-
     G_done_msg(" ");
 
     return 0;

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-07 14:19:08 UTC (rev 44559)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/i.photo.rectify/local_proto.h	2010-12-07 17:49:57 UTC (rev 44560)
@@ -1,5 +1,5 @@
 /* angle.c */
-int camera_angle(struct cache *, char *);
+int camera_angle(char *);
 
 /* aver_z.c */
 int get_aver_elev(struct Ortho_Control_Points *, double *);



More information about the grass-commit mailing list