add alternative bbox forward/inverse functions

AUTHORS

@@ -58,6 +58,7 @@ Patches
  - Philipp Spitzer
  - Dave Stubbs
  - River Tarnell
+ - Oliver Tonnhofer
  - Lennard voor den Dag
  - Shaun Walbridge
  - Leslie Wu

CHANGELOG

@@ -14,6 +14,9 @@ For a complete change history, see the SVN log.
 Mapnik Trunk
 ------------
 
+- Added alternative, more robust proj_transform functions to project a bbox using more points than just the four
+  corners to ensure an optimally sized bbox despite proj4 out of bounds conditions. (olt)
+
 - Added map.base parameter that can be set to control where files with relative paths should be interpreted
   from when a map is loaded from a string or saved to a string. It defaults to an empty string which means
   that the base path will be the current working directory of the mapnik process. When a stylesheet is read

src/proj_transform.cpp

@@ -24,10 +24,15 @@
 
 // mapnik
 #include <mapnik/proj_transform.hpp>
+#include <mapnik/coord.hpp> 
 #include <mapnik/utils.hpp>
+
 // proj4
 #include <proj_api.h>
 
+// stl
+#include <vector>
+
 namespace mapnik {
     
 proj_transform::proj_transform(projection const& source, 
@@ -59,19 +64,19 @@ bool proj_transform::forward (double & x, double & y , double & z) const
 #if defined(MAPNIK_THREADSAFE) && PJ_VERSION < 480
     mutex::scoped_lock lock(projection::mutex_);
 #endif
-        
+
     if (pj_transform( source_.proj_, dest_.proj_, 1, 
                       0, &x,&y,&z) != 0)
     {
         return false;
     }
-        
+
     if (is_dest_longlat_)
     {
         x *= RAD_TO_DEG;
         y *= RAD_TO_DEG;
     }
-        
+
     return true;
 } 
         
@@ -79,13 +84,13 @@ bool proj_transform::backward (double & x, double & y , double & z) const
 {
     if (is_source_equal_dest_)
         return true;
-      
+
     if (is_dest_longlat_)
     {
         x *= DEG_TO_RAD;
         y *= DEG_TO_RAD;
     }
-        
+
 #if defined(MAPNIK_THREADSAFE) && PJ_VERSION < 480
     mutex::scoped_lock lock(projection::mutex_);
 #endif
@@ -95,13 +100,13 @@ bool proj_transform::backward (double & x, double & y , double & z) const
     {
         return false;
     }
-        
+
     if (is_source_longlat_)
     {
         x *= RAD_TO_DEG;
         y *= RAD_TO_DEG;
     }
-        
+
     return true;
 }
 
@@ -109,6 +114,7 @@ bool proj_transform::forward (box2d<double> & box) const
 {
     if (is_source_equal_dest_)
         return true;
+
     double minx = box.minx();
     double miny = box.miny();
     double maxx = box.maxx();
@@ -121,7 +127,7 @@ bool proj_transform::forward (box2d<double> & box) const
     box.init(minx,miny,maxx,maxy);
     return true;
 } 
-        
+
 bool proj_transform::backward (box2d<double> & box) const
 {
     if (is_source_equal_dest_)
@@ -140,6 +146,102 @@ bool proj_transform::backward (box2d<double> & box) const
     return true;
 }
 
+void envelope_points(std::vector< coord<double,2> > & coords, box2d<double>& env, int points)
+{
+    double width = env.width();
+    double height = env.height();
+    
+    int steps;
+    
+    if (points <= 4) {
+        steps = 0;
+    } else {
+        steps = static_cast<int>(ceil((points - 4) / 4.0));
+    }
+    
+    steps += 1;
+    double xstep = width / steps;
+    double ystep = height / steps;
+    
+    for (int i=0; i<=steps; i++) {
+        coords.push_back(coord<double,2>(env.minx() + i * xstep, env.miny()));
+        coords.push_back(coord<double,2>(env.minx() + i * xstep, env.maxy()));
+
+    }
+    for (int i=1; i<steps; i++) {
+        coords.push_back(coord<double,2>(env.minx(), env.miny() + i * ystep));
+        coords.push_back(coord<double,2>(env.maxx(), env.miny() + i * ystep));
+    }
+}
+
+box2d<double> calculate_bbox(std::vector<coord<double,2> > & points) {
+    std::vector<coord<double,2> >::iterator it = points.begin();
+    std::vector<coord<double,2> >::iterator it_end = points.end();
+    
+    box2d<double> env(*it, *(++it));
+    for (; it!=it_end; ++it) {
+        env.expand_to_include(*it);
+    }
+    return env;
+}
+
+
+/* More robust, but expensive, bbox transform
+ * in the face of proj4 out of bounds conditions.
+ * Can result in 20 -> 10 r/s performance hit.
+ * Alternative is to provide proper clipping box
+ * in the target srs by setting map 'maximum-extent' 
+ */
+bool proj_transform::backward(box2d<double>& env, int points) const
+{
+    if (is_source_equal_dest_)
+        return true;
+    
+    std::vector<coord<double,2> > coords;
+    envelope_points(coords, env, points);
+    
+    double z;
+    for (std::vector<coord<double,2> >::iterator it = coords.begin(); it!=coords.end(); ++it) {
+        z = 0;
+        if (!backward(it->x, it->y, z)) {
+            return false;
+        }
+    }
+
+    box2d<double> result = calculate_bbox(coords);
+    
+    env.re_center(result.center().x, result.center().y);
+    env.height(result.height());
+    env.width(result.width());
+
+    return true;
+}
+
+bool proj_transform::forward(box2d<double>& env, int points) const
+{
+    if (is_source_equal_dest_)
+        return true;
+    
+    std::vector<coord<double,2> > coords;
+    envelope_points(coords, env, points);
+    
+    double z;
+    for (std::vector<coord<double,2> >::iterator it = coords.begin(); it!=coords.end(); ++it) {
+        z = 0;
+        if (!forward(it->x, it->y, z)) {
+            return false;
+        }
+    }
+
+    box2d<double> result = calculate_bbox(coords);
+    
+    env.re_center(result.center().x, result.center().y);
+    env.height(result.height());
+    env.width(result.width());
+
+    return true;
+}
+
 mapnik::projection const& proj_transform::source() const
 {
     return source_;