Refactor bordwalk.c

Improve accuracy of -p/-g output


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

raster/r.proj/bordwalk.c

@@ -31,182 +31,203 @@
 #include <grass/raster.h>
 #include <grass/gprojects.h>
 #include <grass/glocale.h>
+#include "r.proj.h"
 
-void bordwalk(struct Cell_head *from_hd, struct Cell_head *to_hd,
-	      struct pj_info *from_pj, struct pj_info *to_pj)
+#ifdef MIN
+#undef MIN
+#endif
+#define MIN(a, b) (((a) < (b)) ? (a) : (b))
+
+#ifdef MAX
+#undef MAX
+#endif
+#define MAX(a, b) (((a) > (b)) ? (a) : (b))
+
+static void debug(const char *name, const struct Cell_head *hd)
 {
-    double idx;
-    double hx, hy;
-    double xmin, xmax;
-    double ymin, ymax;
+    G_debug(3, "%s: xmin: %f; xmax: %f; ymin: %f; ymax: %f",
+	    name, hd->west, hd->east, hd->south, hd->north);
+}
 
-    /* Set some (un)reasonable defaults before we walk the borders */
+static void update(struct Cell_head *to_hd, double hx, double hy)
+{
+    to_hd->east  = MAX(to_hd->east,  hx);
+    to_hd->west  = MIN(to_hd->west,  hx);
+    to_hd->north = MAX(to_hd->north, hy);
+    to_hd->south = MIN(to_hd->south, hy);
+}
 
-    xmax = to_hd->west - 0.000001;
-    xmin = to_hd->east + 0.000001;
-    ymin = to_hd->north + 0.000001;
-    ymax = to_hd->south - 0.000001;
+static void intersect(struct Cell_head *to_hd, const struct Cell_head *from_hd)
+{
+    to_hd->east  = MIN(to_hd->east,  from_hd->east);
+    to_hd->west  = MAX(to_hd->west,  from_hd->west);
+    to_hd->north = MIN(to_hd->north, from_hd->north);
+    to_hd->south = MAX(to_hd->south, from_hd->south);
+}
 
-    /* Start walking */
+static int inside(const struct Cell_head *ref_hd, double hx, double hy)
+{
+    return 
+	(hx <= ref_hd->east) &&
+	(hx >= ref_hd->west) &&
+	(hy <= ref_hd->north) &&
+	(hy >= ref_hd->south);
+}
+
+static void invert(struct Cell_head *cur_hd, const struct Cell_head *ref_hd,
+		   double epsilon)
+{
+    cur_hd->east  = ref_hd->west  - epsilon;
+    cur_hd->west  = ref_hd->east  + epsilon;
+    cur_hd->north = ref_hd->south - epsilon;
+    cur_hd->south = ref_hd->north + epsilon;
+}
+
+static void proj_update(const struct pj_info *from_pj, const struct pj_info *to_pj,
+			struct Cell_head *to_hd, double hx, double hy)
+{
+	if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
+	    return;
+	update(to_hd, hx, hy);
+}
+
+void bordwalk1(const struct pj_info *from_pj, const struct pj_info *to_pj,
+	       const struct Cell_head *from_hd, struct Cell_head *to_hd)
+{
+    double idx;
 
     /* Top */
     for (idx = from_hd->west + from_hd->ew_res / 2; idx < from_hd->east;
-	 idx += from_hd->ew_res) {
-	hx = idx;
-	hy = from_hd->north - from_hd->ns_res / 2;
-	if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
-	    continue;
-	/* check if we are within the region, but allow for some 'almost inside' points */
-	/* (should probably be a factor based on input and output resolutions) */
-	if (!(hx < to_hd->west - to_hd->ew_res) &&
-	    !(hx > to_hd->east + to_hd->ew_res) &&
-	    !(hy < to_hd->south - to_hd->ns_res) &&
-	    !(hy > to_hd->north + to_hd->ns_res)) {
-	    xmin = !(hx > xmin) ? hx : xmin;
-	    xmax = !(hx < xmax) ? hx : xmax;
-	    ymin = !(hy > ymin) ? hy : ymin;
-	    ymax = !(hy < ymax) ? hy : ymax;
-	}
-    }
+	 idx += from_hd->ew_res)
+	proj_update(from_pj, to_pj, to_hd, idx, from_hd->north - from_hd->ns_res / 2);
 
-    G_debug(3, "Top: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin, xmax,
-	    ymin, ymax);
+    debug("Top", to_hd);
 
     /* Right */
     for (idx = from_hd->north - from_hd->ns_res / 2; idx > from_hd->south;
-	 idx -= from_hd->ns_res) {
-	hx = from_hd->east - from_hd->ew_res / 2;
-	hy = idx;
-	if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
-	    continue;
-	if (!(hx < to_hd->west - to_hd->ew_res) &&
-	    !(hx > to_hd->east + to_hd->ew_res) &&
-	    !(hy < to_hd->south - to_hd->ns_res) &&
-	    !(hy > to_hd->north + to_hd->ns_res)) {
-	    xmin = !(hx > xmin) ? hx : xmin;
-	    xmax = !(hx < xmax) ? hx : xmax;
-	    ymin = !(hy > ymin) ? hy : ymin;
-	    ymax = !(hy < ymax) ? hy : ymax;
-	}
-    }
+	 idx -= from_hd->ns_res)
+	proj_update(from_pj, to_pj, to_hd, from_hd->east - from_hd->ew_res / 2, idx);
 
-    G_debug(3, "Right: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin, xmax,
-	    ymin, ymax);
+    debug("Right", to_hd);
 
     /* Bottom */
     for (idx = from_hd->east - from_hd->ew_res / 2; idx > from_hd->west;
-	 idx -= from_hd->ew_res) {
-	hx = idx;
-	hy = from_hd->south + from_hd->ns_res / 2;
-	if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
-	    continue;
-	if (!(hx < to_hd->west - to_hd->ew_res) &&
-	    !(hx > to_hd->east + to_hd->ew_res) &&
-	    !(hy < to_hd->south - to_hd->ns_res) &&
-	    !(hy > to_hd->north + to_hd->ns_res)) {
-	    xmin = !(hx > xmin) ? hx : xmin;
-	    xmax = !(hx < xmax) ? hx : xmax;
-	    ymin = !(hy > ymin) ? hy : ymin;
-	    ymax = !(hy < ymax) ? hy : ymax;
-	}
-    }
+	 idx -= from_hd->ew_res)
+	proj_update(from_pj, to_pj, to_hd, idx, from_hd->south + from_hd->ns_res / 2);
 
-    G_debug(3, "Bottom: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin, xmax,
-	    ymin, ymax);
+    debug("Bottom", to_hd);
 
     /* Left */
     for (idx = from_hd->south + from_hd->ns_res / 2; idx < from_hd->north;
-	 idx += from_hd->ns_res) {
-	hx = from_hd->west + from_hd->ew_res / 2;
-	hy = idx;
-	if (pj_do_proj(&hx, &hy, from_pj, to_pj) < 0)
-	    continue;
-	if (!(hx < to_hd->west - to_hd->ew_res) &&
-	    !(hx > to_hd->east + to_hd->ew_res) &&
-	    !(hy < to_hd->south - to_hd->ns_res) &&
-	    !(hy > to_hd->north + to_hd->ns_res)) {
-	    xmin = !(hx > xmin) ? hx : xmin;
-	    xmax = !(hx < xmax) ? hx : xmax;
-	    ymin = !(hy > ymin) ? hy : ymin;
-	    ymax = !(hy < ymax) ? hy : ymax;
-	}
-    }
+	 idx += from_hd->ns_res)
+	proj_update(from_pj, to_pj, to_hd, from_hd->west + from_hd->ew_res / 2, idx);
 
-    G_debug(3, "Left: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin, xmax,
-	    ymin, ymax);
+    debug("Left", to_hd);
+}
 
-    /* check some special cases by reversing the projection */
+static int proj_inside(const struct pj_info *from_pj, const struct pj_info *to_pj,
+		       const struct Cell_head *ref_hd, double hx, double hy)
+{
+    if (pj_do_proj(&hx, &hy, to_pj, from_pj) < 0)
+	return 0;
+    return inside(ref_hd, hx, hy);
+}
 
-    if (xmin > to_hd->west) {
-	hx = to_hd->west + to_hd->ew_res / 2;
-	hy = to_hd->south + (to_hd->north - to_hd->south) / 2;
-	if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
-	    !(hx < from_hd->west) && !(hx > from_hd->east) &&
-	    !(hy < from_hd->south) && !(hy > from_hd->north))
-	    xmin = to_hd->west + to_hd->ew_res / 2;
+static void reverse_check(const struct pj_info *from_pj, const struct pj_info *to_pj,
+			  const struct Cell_head *from_hd,
+			  const struct Cell_head *to_hd,
+			  struct Cell_head *cur_hd)
+{
+    if (cur_hd->west > to_hd->west) {
+	double hx = to_hd->west + to_hd->ew_res / 2;
+	double hy = to_hd->south + (to_hd->north - to_hd->south) / 2;
+	if (proj_inside(from_pj, to_pj, from_hd, hx, hy))
+	    cur_hd->west = hx;
     }
 
-    if (xmax < to_hd->east) {
-	hx = to_hd->east - to_hd->ew_res / 2;
-	hy = to_hd->south + (to_hd->north - to_hd->south) / 2;
-	if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
-	    !(hx < from_hd->west) && !(hx > from_hd->east) &&
-	    !(hy < from_hd->south) && !(hy > from_hd->north))
-	    xmax = to_hd->east - to_hd->ew_res / 2;
+    if (cur_hd->east < to_hd->east) {
+	double hx = to_hd->east - to_hd->ew_res / 2;
+	double hy = to_hd->south + (to_hd->north - to_hd->south) / 2;
+	if (proj_inside(from_pj, to_pj, from_hd, hx, hy))
+	    cur_hd->east = hx;
     }
 
-    if (ymin > to_hd->south) {
-	hx = to_hd->west + (to_hd->east - to_hd->west) / 2;
-	hy = to_hd->south + to_hd->ns_res / 2;
-	if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
-	    !(hx < from_hd->west) && !(hx > from_hd->east) &&
-	    !(hy < from_hd->south) && !(hy > from_hd->north))
-	    ymin = to_hd->south + to_hd->ns_res / 2;
+    if (cur_hd->south > to_hd->south) {
+	double hx = to_hd->west + (to_hd->east - to_hd->west) / 2;
+	double hy = to_hd->south + to_hd->ns_res / 2;
+	if (proj_inside(from_pj, to_pj, from_hd, hx, hy))
+	    cur_hd->south = hy;
     }
 
-    if (ymax < to_hd->north) {
-	hx = to_hd->west + (to_hd->east - to_hd->west) / 2;
-	hy = to_hd->north - to_hd->ns_res / 2;
-	if (!(pj_do_proj(&hx, &hy, to_pj, from_pj) < 0) &&
-	    !(hx < from_hd->west) && !(hx > from_hd->east) &&
-	    !(hy < from_hd->south) && !(hy > from_hd->north))
-	    ymax = to_hd->north - to_hd->ns_res / 2;
+    if (cur_hd->north < to_hd->north) {
+	double hx = to_hd->west + (to_hd->east - to_hd->west) / 2;
+	double hy = to_hd->north - to_hd->ns_res / 2;
+	if (proj_inside(from_pj, to_pj, from_hd, hx, hy))
+	    cur_hd->north = hy;
     }
+}
 
-    G_debug(3, "Extra check: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin,
-	    xmax, ymin, ymax);
+static int outside(const struct Cell_head *cur_hd, const struct Cell_head *ref_hd)
+{
+    return
+	(cur_hd->west  > ref_hd->east ) ||
+	(cur_hd->east  < ref_hd->west ) ||
+	(cur_hd->south > ref_hd->north) ||
+	(cur_hd->north < ref_hd->south);
+}
+
+static void snap_to_grid(struct Cell_head *cur_hd, const struct Cell_head *ref_hd)
+{
+    int lidx = (int) floor(Rast_easting_to_col( cur_hd->west,  ref_hd));
+    int ridx = (int) floor(Rast_easting_to_col( cur_hd->east,  ref_hd));
+    int bidx = (int) floor(Rast_northing_to_row(cur_hd->south, ref_hd));
+    int tidx = (int) floor(Rast_northing_to_row(cur_hd->north, ref_hd));
+
+    cur_hd->west  = Rast_col_to_easting( lidx + 0.0, ref_hd);
+    cur_hd->east  = Rast_col_to_easting( ridx + 1.0, ref_hd);
+    cur_hd->south = Rast_row_to_northing(bidx + 1.0, ref_hd);
+    cur_hd->north = Rast_row_to_northing(tidx + 0.0, ref_hd);
+}
+
+void bordwalk(const struct Cell_head *from_hd, struct Cell_head *to_hd,
+	      const struct pj_info *from_pj, const struct pj_info *to_pj)
+{
+    struct Cell_head cur_hd;
+
+    /* Set some (un)reasonable defaults before we walk the borders */
+
+    invert(&cur_hd, to_hd, 1.0e-6);
+
+    /* Start walking */
+
+    bordwalk1(from_pj, to_pj, from_hd, &cur_hd);
+
+    intersect(&cur_hd, to_hd);
+
+    /* check some special cases by reversing the projection */
+
+    reverse_check(from_pj, to_pj, from_hd, to_hd, &cur_hd);
+
+    debug("Extra check", &cur_hd);
 
     /* if we still have some unresonable default minmax left, then abort */
 
-    if ((xmin > to_hd->east) || (xmax < to_hd->west)
-	|| (ymin > to_hd->north) || (ymax < to_hd->south))
+    if (outside(&cur_hd, to_hd))
 	G_fatal_error(_("Input raster map is outside current region"));
 
-    if (xmin < to_hd->west + to_hd->ew_res / 2)
-	xmin = to_hd->west + to_hd->ew_res / 2;
-    if (xmax > to_hd->east - to_hd->ew_res / 2)
-	xmax = to_hd->east - to_hd->ew_res / 2;
-    if (ymin < to_hd->south + to_hd->ns_res / 2)
-	ymin = to_hd->south + to_hd->ns_res / 2;
-    if (ymax > to_hd->north - to_hd->ns_res / 2)
-	ymax = to_hd->north - to_hd->ns_res / 2;
+    intersect(&cur_hd, to_hd);
 
     /* adjust to edges */
 
-    idx = (int)floor(Rast_easting_to_col(xmin, to_hd));
-    xmin = Rast_col_to_easting(idx + 0.0, to_hd);
-    idx = (int)floor(Rast_easting_to_col(xmax, to_hd));
-    xmax = Rast_col_to_easting(idx + 1.0, to_hd);
-    idx = (int)floor(Rast_northing_to_row(ymin, to_hd));
-    ymin = Rast_row_to_northing(idx + 1.0, to_hd);
-    idx = (int)floor(Rast_northing_to_row(ymax, to_hd));
-    ymax = Rast_row_to_northing(idx + 0.0, to_hd);
-
-    to_hd->west = (xmin < to_hd->west) ? to_hd->west : xmin;
-    to_hd->east = (xmax > to_hd->east) ? to_hd->east : xmax;
-    to_hd->south = (ymin < to_hd->south) ? to_hd->south : ymin;
-    to_hd->north = (ymax > to_hd->north) ? to_hd->north : ymax;
-
-    G_debug(3, "Final check: xmin: %f; xmax: %f; ymin: %f; ymax: %f", xmin,
-	    xmax, ymin, ymax);
+    snap_to_grid(&cur_hd, to_hd);
+
+    intersect(to_hd, &cur_hd);
+
+    debug("Final check", to_hd);
+}
+
+void bordwalk2(const struct Cell_head *from_hd, struct Cell_head *to_hd,
+	      const struct pj_info *from_pj, const struct pj_info *to_pj)
+{
+    bordwalk1(from_pj, to_pj, from_hd, to_hd);
 }

raster/r.proj/main.c

@@ -245,6 +245,7 @@ int main(int argc, char **argv)
 
     mapname = outmap->answer ? outmap->answer : inmap->answer;
     if (mapname && !list->answer && !overwrite &&
+	!print_bounds->answer && !gprint_bounds->answer &&
 	G_find_raster(mapname, G_mapset()))
 	G_fatal_error(_("option <%s>: <%s> exists."), "output", mapname);
 
@@ -354,10 +355,15 @@ int main(int argc, char **argv)
 	G_message(_("Input map <%s@%s> in location <%s>:"),
 	    inmap->answer, setname, inlocation->answer);
 
-	if (pj_do_proj(&iwest, &isouth, &iproj, &oproj) < 0)
-	    G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
-	if (pj_do_proj(&ieast, &inorth, &iproj, &oproj) < 0)
-	    G_fatal_error(_("Error in pj_do_proj (projection of input coordinate pair)"));
+	outcellhd.north = -1e9;
+	outcellhd.south =  1e9;
+	outcellhd.east  = -1e9;
+	outcellhd.west  =  1e9;
+	bordwalk2(&incellhd, &outcellhd, &iproj, &oproj);
+	inorth = outcellhd.north;
+	isouth = outcellhd.south;
+	ieast  = outcellhd.east;
+	iwest  = outcellhd.west;
 
 	G_format_northing(inorth, north_str, curr_proj);
 	G_format_northing(isouth, south_str, curr_proj);

raster/r.proj/r.proj.h

@@ -35,8 +35,10 @@ struct menu
     char *text;			/* menu display - full description       */
 };
 
-extern void bordwalk(struct Cell_head *, struct Cell_head *, struct pj_info *,
-		     struct pj_info *);
+extern void bordwalk(const struct Cell_head *, struct Cell_head *,
+		     const struct pj_info *, const struct pj_info *);
+extern void bordwalk2(const struct Cell_head *, struct Cell_head *,
+		      const struct pj_info *, const struct pj_info *);
 extern struct cache *readcell(int, const char *);
 extern block *get_block(struct cache *, int);
 extern void release_cache(struct cache *);