simgear/bucket/newbucket.cxx

   1 /**************************************************************************
   2  * newbucket.hxx -- new bucket routines for better world modeling
   3  *
   4  * Written by Curtis L. Olson, started February 1999.
   5  *
   6  * Copyright (C) 1999  Curtis L. Olson - http://www.flightgear.org/~curt/
   7  *
   8  * This library is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU Library General Public
  10  * License as published by the Free Software Foundation; either
  11  * version 2 of the License, or (at your option) any later version.
  12  *
  13  * This library is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16  * Library General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU General Public License
  19  * along with this program; if not, write to the Free Software
  20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  21  *
  22  * $Id$
  23  **************************************************************************/
  24
  25
  26 #ifdef HAVE_CONFIG_H
  27 #  include <simgear_config.h>
  28 #endif
  29
  30 #include <cmath>
  31 #include <cstdio> // some platforms need this for ::snprintf
  32 #include <iostream>
  33
  34 #include <simgear/misc/sg_path.hxx>
  35 #include <simgear/debug/logstream.hxx>
  36
  37 #include "newbucket.hxx"
  38
  39
  40 // default constructor
  41 SGBucket::SGBucket() :
  42     lon(-1000),
  43     lat(-1000),
  44     x(0),
  45     y(0)
  46 {
  47 }
  48
  49 bool SGBucket::isValid() const
  50 {
  51     // The most northerly valid latitude is 89, not 90. There is no tile
  52     // whose *bottom* latitude is 90. Similar there is no tile whose left egde
  53     // is 180 longitude.
  54     return (lon >= -180) &&
  55             (lon < 180) &&
  56             (lat >= -90) &&
  57             (lat < 90) &&
  58             (x < 8) && (y < 8);
  59 }
  60
  61 void SGBucket::make_bad()
  62 {
  63     lon = -1000;
  64     lat = -1000;
  65 }
  66
  67 // constructor for specified location
  68 SGBucket::SGBucket(const double dlon, const double dlat) {
  69     set_bucket(dlon, dlat);
  70 }
  71
  72 SGBucket::SGBucket(const SGGeod& geod) {
  73     set_bucket(geod.getLongitudeDeg(),
  74                    geod.getLatitudeDeg());
  75 }
  76
  77 // Parse a unique scenery tile index and find the lon, lat, x, and y
  78 SGBucket::SGBucket(const long int bindex) {
  79     long int index = bindex;
  80
  81     lon = index >> 14;
  82     index -= lon << 14;
  83     lon -= 180;
  84
  85     lat = index >> 6;
  86     index -= lat << 6;
  87     lat -= 90;
  88
  89     y = index >> 3;
  90     index -= y << 3;
  91
  92     x = index;
  93 }
  94
  95 /* Calculate the greatest integral value less than
  96  * or equal to the given value (floor(x)),
  97  * but attribute coordinates close to the boundary to the next
  98  * (increasing) integral
  99  */
 100 static int floorWithEpsilon(double x)
 101 {
 102     double diff = x - static_cast<int>(x);
 103     if ( (x >= 0.0) || (fabs(diff) < SG_EPSILON) ) {
 104         return static_cast<int>(x);
 105     } else {
 106         return static_cast<int>(x) - 1;
 107     }
 108 }
 109
 110 // Set the bucket params for the specified lat and lon
 111 void SGBucket::set_bucket( double dlon, double dlat )
 112 {
 113     if ((dlon < -180.0) || (dlon >= 180.0)) {
 114         SG_LOG(SG_TERRAIN, SG_WARN, "SGBucket::set_bucket: passed longitude:" << dlon);
 115         dlon = SGMiscd::normalizePeriodic(-180.0, 180.0, dlon);
 116     }
 117
 118     if ((dlat < -90.0) || (dlat > 90.0)) {
 119         SG_LOG(SG_TERRAIN, SG_WARN, "SGBucket::set_bucket: passed latitude" << dlat);
 120         dlat = SGMiscd::clip(dlat, -90.0, 90.0);
 121     }
 122
 123     //
 124     // longitude first
 125     //
 126     double span = sg_bucket_span( dlat );
 127     // we do NOT need to special case lon=180 here, since
 128     // normalizePeriodic will never return 180; it will
 129     // return -180, which is what we want.
 130     lon = floorWithEpsilon(dlon);
 131
 132     // find subdivision or super lon if needed
 133     if ( span <= 1.0 ) {
 134         /* We have more than one tile per degree of
 135          * longitude, so we need an x offset.
 136          */
 137         x = (int)((dlon - lon) / span);
 138     } else {
 139         /* We have one or more degrees per tile,
 140          * so we need to find the base longitude
 141          * of that tile.
 142          *
 143          * First we calculate the integral base longitude
 144          * (e.g. -85.5 => -86) and then find the greatest
 145          * multiple of span that is less than or equal to
 146          * that longitude.
 147          *
 148          * That way, the Greenwich Meridian is always
 149          * a tile border.
 150          */
 151         lon=static_cast<int>(floor(lon / span) * span);
 152         x = 0;
 153     }
 154
 155     //
 156     // then latitude
 157     //
 158     lat = floorWithEpsilon(dlat);
 159
 160     // special case when passing in the north pole point (possibly due to
 161     // clipping latitude above). Ensures we generate a valid bucket in this
 162     // scenario
 163     if (lat == 90) {
 164         lat = 89;
 165         y = 7;
 166     } else {
 167         /* Latitude base and offset are easier, as
 168          * tiles always are 1/8 degree of latitude wide.
 169          */
 170         y = (int)((dlat - lat) * 8);
 171     }
 172 }
 173
 174 void SGBucket::set_bucket(const SGGeod& geod)
 175 {
 176     set_bucket(geod.getLongitudeDeg(), geod.getLatitudeDeg());
 177 }
 178
 179 // Build the path name for this bucket
 180 std::string SGBucket::gen_base_path() const {
 181     // long int index;
 182     int top_lon, top_lat, main_lon, main_lat;
 183     char hem, pole;
 184     char raw_path[256];
 185
 186     top_lon = lon / 10;
 187     main_lon = lon;
 188     if ( (lon < 0) && (top_lon * 10 != lon) ) {
 189         top_lon -= 1;
 190     }
 191     top_lon *= 10;
 192     if ( top_lon >= 0 ) {
 193         hem = 'e';
 194     } else {
 195         hem = 'w';
 196         top_lon *= -1;
 197     }
 198     if ( main_lon < 0 ) {
 199         main_lon *= -1;
 200     }
 201
 202     top_lat = lat / 10;
 203     main_lat = lat;
 204     if ( (lat < 0) && (top_lat * 10 != lat) ) {
 205         top_lat -= 1;
 206     }
 207     top_lat *= 10;
 208     if ( top_lat >= 0 ) {
 209         pole = 'n';
 210     } else {
 211         pole = 's';
 212         top_lat *= -1;
 213     }
 214     if ( main_lat < 0 ) {
 215         main_lat *= -1;
 216     }
 217
 218     ::snprintf(raw_path, 256, "%c%03d%c%02d/%c%03d%c%02d",
 219             hem, top_lon, pole, top_lat,
 220             hem, main_lon, pole, main_lat);
 221
 222     SGPath path( raw_path );
 223
 224     return path.str();
 225 }
 226
 227
 228 // return width of the tile in degrees
 229 double SGBucket::get_width() const {
 230     return sg_bucket_span( get_center_lat() );
 231 }
 232
 233
 234 // return height of the tile in degrees
 235 double SGBucket::get_height() const {
 236     return SG_BUCKET_SPAN;
 237 }
 238
 239 double SGBucket::get_highest_lat() const
 240 {
 241     unsigned char adjustedY = y;
 242     if (lat >= 0) {
 243         // tile is north of the equator, so we want the top edge. Add one
 244         // to y to achieve this.
 245         ++adjustedY;
 246     }
 247
 248         return lat + (adjustedY / 8.0);
 249 }
 250
 251
 252 // return width of the tile in meters. This function is used by the
 253 // tile-manager to estimate how many tiles are in the view distance, so
 254 // we care about the smallest width, which occurs at the highest latitude.
 255 double SGBucket::get_width_m() const
 256 {
 257     double clat_rad = get_highest_lat() * SGD_DEGREES_TO_RADIANS;
 258     double cos_lat = cos( clat_rad );
 259     if (fabs(cos_lat) < SG_EPSILON) {
 260         // happens for polar tiles, since we pass in a latitude of 90
 261         // return an arbitrary small value so all tiles are loaded
 262         return 10.0;
 263     }
 264
 265     double local_radius = cos_lat * SG_EQUATORIAL_RADIUS_M;
 266     double local_perimeter = local_radius * SGD_2PI;
 267     double degree_width = local_perimeter / 360.0;
 268
 269     return get_width() * degree_width;
 270 }
 271
 272
 273 // return height of the tile in meters
 274 double SGBucket::get_height_m() const {
 275     double perimeter = SG_EQUATORIAL_RADIUS_M * SGD_2PI;
 276     double degree_height = perimeter / 360.0;
 277
 278     return SG_BUCKET_SPAN * degree_height;
 279 }
 280
 281 SGBucket SGBucket::sibling(int dx, int dy) const
 282 {
 283     if (!isValid()) {
 284         SG_LOG(SG_TERRAIN, SG_WARN, "SGBucket::sibling: requesting sibling of invalid bucket");
 285         return SGBucket();
 286     }
 287
 288     double clat = get_center_lat() + dy * SG_BUCKET_SPAN;
 289     // return invalid here instead of clipping, so callers can discard
 290     // invalid buckets without having to check if it's an existing one
 291     if ((clat < -90.0) || (clat > 90.0)) {
 292         return SGBucket();
 293     }
 294
 295     // find the lon span for the new latitude
 296     double span = sg_bucket_span( clat );
 297
 298     double tmp = get_center_lon() + dx * span;
 299     tmp = SGMiscd::normalizePeriodic(-180.0, 180.0, tmp);
 300     return SGBucket(tmp, clat);
 301 }
 302
 303 std::string SGBucket::gen_index_str() const
 304 {
 305         char tmp[20];
 306         ::snprintf(tmp, 20, "%ld",
 307                  (((long)lon + 180) << 14) + ((lat + 90) << 6)
 308                  + (y << 3) + x);
 309         return (std::string)tmp;
 310 }
 311
 312 // find the bucket which is offset by the specified tile units in the
 313 // X & Y direction.  We need the current lon and lat to resolve
 314 // ambiguities when going from a wider tile to a narrower one above or
 315 // below.  This assumes that we are feeding in
 316 SGBucket sgBucketOffset( double dlon, double dlat, int dx, int dy ) {
 317     SGBucket result( dlon, dlat );
 318     double clat = result.get_center_lat() + dy * SG_BUCKET_SPAN;
 319
 320     // walk dy units in the lat direction
 321     result.set_bucket( dlon, clat );
 322
 323     // find the lon span for the new latitude
 324     double span = sg_bucket_span( clat );
 325
 326     // walk dx units in the lon direction
 327     double tmp = dlon + dx * span;
 328     while ( tmp < -180.0 ) {
 329         tmp += 360.0;
 330     }
 331     while ( tmp >= 180.0 ) {
 332         tmp -= 360.0;
 333     }
 334     result.set_bucket( tmp, clat );
 335
 336     return result;
 337 }
 338
 339
 340 // calculate the offset between two buckets
 341 void sgBucketDiff( const SGBucket& b1, const SGBucket& b2, int *dx, int *dy ) {
 342
 343     // Latitude difference
 344     double c1_lat = b1.get_center_lat();
 345     double c2_lat = b2.get_center_lat();
 346     double diff_lat = c2_lat - c1_lat;
 347
 348 #ifdef HAVE_RINT
 349     *dy = (int)rint( diff_lat / SG_BUCKET_SPAN );
 350 #else
 351     if ( diff_lat > 0 ) {
 352         *dy = (int)( diff_lat / SG_BUCKET_SPAN + 0.5 );
 353     } else {
 354         *dy = (int)( diff_lat / SG_BUCKET_SPAN - 0.5 );
 355     }
 356 #endif
 357
 358     // longitude difference
 359     double diff_lon=0.0;
 360     double span=0.0;
 361
 362     SGBucket tmp_bucket;
 363     // To handle crossing the bucket size boundary
 364     //  we need to account for different size buckets.
 365
 366     if ( sg_bucket_span(c1_lat) <= sg_bucket_span(c2_lat) )
 367     {
 368         span = sg_bucket_span(c1_lat);
 369     } else {
 370         span = sg_bucket_span(c2_lat);
 371     }
 372
 373     diff_lon = b2.get_center_lon() - b1.get_center_lon();
 374
 375     if (diff_lon <0.0)
 376     {
 377        diff_lon -= b1.get_width()*0.5 + b2.get_width()*0.5 - span;
 378     }
 379     else
 380     {
 381        diff_lon += b1.get_width()*0.5 + b2.get_width()*0.5 - span;
 382     }
 383
 384
 385 #ifdef HAVE_RINT
 386     *dx = (int)rint( diff_lon / span );
 387 #else
 388     if ( diff_lon > 0 ) {
 389         *dx = (int)( diff_lon / span + 0.5 );
 390     } else {
 391         *dx = (int)( diff_lon / span - 0.5 );
 392     }
 393 #endif
 394 }
 395
 396 void sgGetBuckets( const SGGeod& min, const SGGeod& max, std::vector<SGBucket>& list ) {
 397     double lon, lat, span;
 398
 399     for (lat = min.getLatitudeDeg(); lat <= max.getLatitudeDeg(); lat += SG_BUCKET_SPAN) {
 400         span = sg_bucket_span( lat );
 401         for (lon = min.getLongitudeDeg(); lon <= max.getLongitudeDeg(); lon += span)
 402         {
 403             SGBucket b(lon, lat);
 404             if (!b.isValid()) {
 405                 continue;
 406             }
 407
 408             list.push_back(b);
 409         }
 410     }
 411 }
 412
 413 std::ostream& operator<< ( std::ostream& out, const SGBucket& b )
 414 {
 415     return out << b.lon << ":" << (int)b.x << ", " << b.lat << ":" << (int)b.y;
 416 }
 417