add resolution option, more robust resolution calculation

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@43886 15284696-431f-4ddb-bdfa-cd5b030d7da7

imagery/i.rectify/get_wind.c

@@ -1,16 +1,25 @@
+#include <math.h>
 #include <grass/glocale.h>
 #include "global.h"
 #include "crs.h"		/* CRS HEADER FILE */
 
-int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order)
+int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order, double res)
 {
-    double n, e;
+    double n, e, ad;
+    struct _corner {
+        double n, e;
+    } nw, ne, se, sw;
 
+    /* extends */
     CRS_georef(w1->west, w1->north, &e, &n, E12, N12, order);
     w2->north = w2->south = n;
     w2->west = w2->east = e;
+    nw.n = n;
+    nw.e = e;
 
     CRS_georef(w1->east, w1->north, &e, &n, E12, N12, order);
+    ne.n = n;
+    ne.e = e;
     if (n > w2->north)
 	w2->north = n;
     if (n < w2->south)
@@ -21,6 +30,8 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order)
 	w2->west = e;
 
     CRS_georef(w1->west, w1->south, &e, &n, E12, N12, order);
+    sw.n = n;
+    sw.e = e;
     if (n > w2->north)
 	w2->north = n;
     if (n < w2->south)
@@ -31,6 +42,8 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order)
 	w2->west = e;
 
     CRS_georef(w1->east, w1->south, &e, &n, E12, N12, order);
+    se.n = n;
+    se.e = e;
     if (n > w2->north)
 	w2->north = n;
     if (n < w2->south)
@@ -40,9 +53,63 @@ int georef_window(struct Cell_head *w1, struct Cell_head *w2, int order)
     if (e < w2->west)
 	w2->west = e;
 
-    w2->ns_res = (w2->north - w2->south) / w1->rows;
-    w2->ew_res = (w2->east - w2->west) / w1->cols;
+    /* resolution */
+    if (res > 0)
+	w2->ew_res = w2->ns_res = res;
+    else {
+	/* this results in ugly res values, and ns_res != ew_res */
+	/* and is no good for rotation */
+	/*
+	w2->ns_res = (w2->north - w2->south) / w1->rows;
+	w2->ew_res = (w2->east - w2->west) / w1->cols;
+	*/
 
+	/* alternative: account for rotation and order > 1 */
+
+	/* N-S extends along western and eastern edge */
+	w2->ns_res = (sqrt((nw.n - sw.n) * (nw.n - sw.n) +
+			  (nw.e - sw.e) * (nw.e - sw.e)) +
+		     sqrt((ne.n - se.n) * (ne.n - se.n) +
+			  (ne.e - se.e) * (ne.e - se.e))) / (2.0 * w1->rows);
+
+	/* E-W extends along northern and southern edge */
+	w2->ew_res = (sqrt((nw.n - ne.n) * (nw.n - ne.n) +
+			  (nw.e - ne.e) * (nw.e - ne.e)) +
+		     sqrt((sw.n - se.n) * (sw.n - se.n) +
+			  (sw.e - se.e) * (sw.e - se.e))) / (2.0 * w1->cols);
+
+	/* make ew_res = ns_res */
+	w2->ns_res = (w2->ns_res + w2->ew_res) / 2.0;
+	w2->ew_res = w2->ns_res;
+
+	/* nice round values */
+	if (w2->ns_res > 1) {
+	    if (w2->ns_res < 10) {
+		/* round to first decimal */
+		w2->ns_res = (int)(w2->ns_res * 10 + 0.5) / 10.0;
+		w2->ew_res = w2->ns_res;
+	    }
+	    else {
+		/* round to integer */
+		w2->ns_res = (int)(w2->ns_res + 0.5);
+		w2->ew_res = w2->ns_res;
+	    }
+	}
+    }
+
+    /* adjust extends */
+    ad = w2->north > 0 ? 0.5 : -0.5;
+    w2->north = (int) (ceil(w2->north / w2->ns_res) + ad) * w2->ns_res;
+    ad = w2->south > 0 ? 0.5 : -0.5;
+    w2->south = (int) (floor(w2->south / w2->ns_res) + ad) * w2->ns_res;
+    ad = w2->east > 0 ? 0.5 : -0.5;
+    w2->east = (int) (ceil(w2->east / w2->ew_res) + ad) * w2->ew_res;
+    ad = w2->west > 0 ? 0.5 : -0.5;
+    w2->west = (int) (floor(w2->west / w2->ew_res) + ad) * w2->ew_res;
+
+    w2->rows = (w2->north - w2->south + w2->ns_res / 2.0) / w2->ns_res;
+    w2->cols = (w2->east - w2->west + w2->ew_res / 2.0) / w2->ew_res;
+    
     G_message(_("Region N=%f S=%f E=%f W=%f"), w2->north, w2->south,
 	    w2->east, w2->west);
     G_message(_("Resolution EW=%f NS=%f"), w2->ew_res, w2->ns_res);

imagery/i.rectify/global.h

@@ -55,7 +55,7 @@ int exec_rectify(int, char *);
 
 /* get_wind.c */
 int get_target_window(int);
-int georef_window(struct Cell_head *, struct Cell_head *, int);
+int georef_window(struct Cell_head *, struct Cell_head *, int, double);
 
 /* matrix.c */
 int compute_georef_matrix(struct Cell_head *, struct Cell_head *, int);

imagery/i.rectify/main.c

@@ -69,10 +69,10 @@ int main(int argc, char *argv[])
     char group[INAME_LEN], extension[INAME_LEN];
     char result[NFILES][15];
     int order;			/* ADDED WITH CRS MODIFICATIONS */
-    int n, i, m, k;
+    int n, i, m, k = 0;
     int got_file = 0;
 
-    struct Option *grp, *val, *ifile, *ext;
+    struct Option *grp, *val, *ifile, *ext, *tres;
     struct Flag *c, *a;
     struct GModule *module;
 
@@ -109,6 +109,12 @@ int main(int argc, char *argv[])
     val->required = YES;
     val->description = _("Rectification polynom order (1-3)");
 
+    tres = G_define_option();
+    tres->key = "res";
+    tres->type = TYPE_DOUBLE;
+    tres->required = NO;
+    tres->description = _("Target resolution (ignored if -c flag used)");
+
     c = G_define_flag();
     c->key = 'c';
     c->description =
@@ -143,8 +149,11 @@ int main(int argc, char *argv[])
 		      MAXORDER);
 
     /* determine the number of files in this group */
-    if (I_get_group_ref(group, &ref) <= 0)
+    if (I_get_group_ref(group, &ref) <= 0) {
+	G_warning(_("Location: %s"), G_location());
+	G_warning(_("Mapset: %s"), G_mapset());
 	G_fatal_error(_("Group <%s> does not exist"), grp->answer);
+    }
 
     if (ref.nfiles <= 0) {
 	G_important_message(_("Group <%s> contains no raster maps; run i.group"),
@@ -184,21 +193,20 @@ int main(int argc, char *argv[])
     /* get the target */
     get_target(group);
 
-
-    if (c->answer) {
-	/* Use current Region */
-	select_target_env();
-	G__get_window(&target_window, "", "WIND", G_mapset());
-	select_current_env();
-    }
-    else {
+    if (!c->answer) {
+	double res = -1;
+	
+	if (tres->answer) {
+	    if (!((res = atof(tres->answer)) > 0))
+		G_warning(_("Target resolution must be > 0, ignored"));
+	}
 	/* Calculate smallest region */
 	if (a->answer)
 	    Rast_get_cellhd(ref.file[0].name, ref.file[0].mapset, &cellhd);
 	else
 	    Rast_get_cellhd(ifile->answers[0], ref.file[0].mapset, &cellhd);
 
-	georef_window(&cellhd, &target_window, order);
+	georef_window(&cellhd, &target_window, order, res);
     }
 
     G_verbose_message(_("Using region: N=%f S=%f, E=%f W=%f"), target_window.north,

imagery/i.rectify/rowcol.h

@@ -1 +1 @@
-#define ROWCOL short
+#define ROWCOL int