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 #ifndef NO_DEPRECATED_API
  68
  69 // constructor for specified location
  70 SGBucket::SGBucket(const double dlon, const double dlat) {
  71     set_bucket(dlon, dlat);
  72 }
  73 #endif
  74
  75 SGBucket::SGBucket(const SGGeod& geod) {
  76     innerSet(geod.getLongitudeDeg(),
  77              geod.getLatitudeDeg());
  78 }
  79
  80 // Parse a unique scenery tile index and find the lon, lat, x, and y
  81 SGBucket::SGBucket(const long int bindex) {
  82     long int index = bindex;
  83
  84     lon = index >> 14;
  85     index -= lon << 14;
  86     lon -= 180;
  87
  88     lat = index >> 6;
  89     index -= lat << 6;
  90     lat -= 90;
  91
  92     y = index >> 3;
  93     index -= y << 3;
  94
  95     x = index;
  96 }
  97
  98 /* Calculate the greatest integral value less than
  99  * or equal to the given value (floor(x)),
 100  * but attribute coordinates close to the boundary to the next
 101  * (increasing) integral
 102  */
 103 static int floorWithEpsilon(double x)
 104 {
 105     return static_cast<int>(floor(x + SG_EPSILON));
 106 }
 107
 108 #ifndef NO_DEPRECATED_API
 109
 110 void SGBucket::set_bucket(double dlon, double dlat)
 111 {
 112     innerSet(dlon, dlat);
 113 }
 114
 115
 116 void SGBucket::set_bucket(const SGGeod& geod)
 117 {
 118     innerSet(geod.getLongitudeDeg(), geod.getLatitudeDeg());
 119 }
 120
 121 #endif
 122
 123 // Set the bucket params for the specified lat and lon
 124 void SGBucket::innerSet( double dlon, double dlat )
 125 {
 126     if ((dlon < -180.0) || (dlon >= 180.0)) {
 127         SG_LOG(SG_TERRAIN, SG_WARN, "SGBucket::set_bucket: passed longitude:" << dlon);
 128         dlon = SGMiscd::normalizePeriodic(-180.0, 180.0, dlon);
 129     }
 130
 131     if ((dlat < -90.0) || (dlat > 90.0)) {
 132         SG_LOG(SG_TERRAIN, SG_WARN, "SGBucket::set_bucket: passed latitude" << dlat);
 133         dlat = SGMiscd::clip(dlat, -90.0, 90.0);
 134     }
 135
 136     //
 137     // longitude first
 138     //
 139     double span = sg_bucket_span( dlat );
 140     // we do NOT need to special case lon=180 here, since
 141     // normalizePeriodic will never return 180; it will
 142     // return -180, which is what we want.
 143     lon = floorWithEpsilon(dlon);
 144
 145     // find subdivision or super lon if needed
 146     if ( span <= 1.0 ) {
 147         /* We have more than one tile per degree of
 148          * longitude, so we need an x offset.
 149          */
 150         x = floorWithEpsilon((dlon - lon) / span);
 151     } else {
 152         /* We have one or more degrees per tile,
 153          * so we need to find the base longitude
 154          * of that tile.
 155          *
 156          * First we calculate the integral base longitude
 157          * (e.g. -85.5 => -86) and then find the greatest
 158          * multiple of span that is less than or equal to
 159          * that longitude.
 160          *
 161          * That way, the Greenwich Meridian is always
 162          * a tile border.
 163          */
 164         lon=static_cast<int>(floor(lon / span) * span);
 165         x = 0;
 166     }
 167
 168     //
 169     // then latitude
 170     //
 171     lat = floorWithEpsilon(dlat);
 172
 173     // special case when passing in the north pole point (possibly due to
 174     // clipping latitude above). Ensures we generate a valid bucket in this
 175     // scenario
 176     if (lat == 90) {
 177         lat = 89;
 178         y = 7;
 179     } else {
 180         /* Latitude base and offset are easier, as
 181          * tiles always are 1/8 degree of latitude wide.
 182          */
 183         y = floorWithEpsilon((dlat - lat) * 8);
 184     }
 185 }
 186
 187 // Build the path name for this bucket
 188 std::string SGBucket::gen_base_path() const {
 189     // long int index;
 190     int top_lon, top_lat, main_lon, main_lat;
 191     char hem, pole;
 192     char raw_path[256];
 193
 194     top_lon = lon / 10;
 195     main_lon = lon;
 196     if ( (lon < 0) && (top_lon * 10 != lon) ) {
 197         top_lon -= 1;
 198     }
 199     top_lon *= 10;
 200     if ( top_lon >= 0 ) {
 201         hem = 'e';
 202     } else {
 203         hem = 'w';
 204         top_lon *= -1;
 205     }
 206     if ( main_lon < 0 ) {
 207         main_lon *= -1;
 208     }
 209
 210     top_lat = lat / 10;
 211     main_lat = lat;
 212     if ( (lat < 0) && (top_lat * 10 != lat) ) {
 213         top_lat -= 1;
 214     }
 215     top_lat *= 10;
 216     if ( top_lat >= 0 ) {
 217         pole = 'n';
 218     } else {
 219         pole = 's';
 220         top_lat *= -1;
 221     }
 222     if ( main_lat < 0 ) {
 223         main_lat *= -1;
 224     }
 225
 226     ::snprintf(raw_path, 256, "%c%03d%c%02d/%c%03d%c%02d",
 227             hem, top_lon, pole, top_lat,
 228             hem, main_lon, pole, main_lat);
 229
 230     SGPath path( raw_path );
 231
 232     return path.str();
 233 }
 234
 235
 236 // return width of the tile in degrees
 237 double SGBucket::get_width() const {
 238     return sg_bucket_span( get_center_lat() );
 239 }
 240
 241
 242 // return height of the tile in degrees
 243 double SGBucket::get_height() const {
 244     return SG_BUCKET_SPAN;
 245 }
 246
 247 double SGBucket::get_highest_lat() const
 248 {
 249     unsigned char adjustedY = y;
 250     if (lat >= 0) {
 251         // tile is north of the equator, so we want the top edge. Add one
 252         // to y to achieve this.
 253         ++adjustedY;
 254     }
 255
 256         return lat + (adjustedY / 8.0);
 257 }
 258
 259
 260 // return width of the tile in meters. This function is used by the
 261 // tile-manager to estimate how many tiles are in the view distance, so
 262 // we care about the smallest width, which occurs at the highest latitude.
 263 double SGBucket::get_width_m() const
 264 {
 265     double clat_rad = get_highest_lat() * SGD_DEGREES_TO_RADIANS;
 266     double cos_lat = cos( clat_rad );
 267     if (fabs(cos_lat) < SG_EPSILON) {
 268         // happens for polar tiles, since we pass in a latitude of 90
 269         // return an arbitrary small value so all tiles are loaded
 270         return 10.0;
 271     }
 272
 273     double local_radius = cos_lat * SG_EQUATORIAL_RADIUS_M;
 274     double local_perimeter = local_radius * SGD_2PI;
 275     double degree_width = local_perimeter / 360.0;
 276
 277     return get_width() * degree_width;
 278 }
 279
 280
 281 // return height of the tile in meters
 282 double SGBucket::get_height_m() const {
 283     double perimeter = SG_EQUATORIAL_RADIUS_M * SGD_2PI;
 284     double degree_height = perimeter / 360.0;
 285
 286     return SG_BUCKET_SPAN * degree_height;
 287 }
 288
 289 unsigned int SGBucket::siblings( int dx, int dy, std::vector<SGBucket>& buckets ) const
 290 {
 291     if (!isValid()) {
 292         SG_LOG(SG_TERRAIN, SG_WARN, "SGBucket::sibling: requesting sibling of invalid bucket");
 293         return 0;
 294     }
 295
 296     double src_span = sg_bucket_span( get_center_lat() );
 297
 298     double clat = get_center_lat() + dy * SG_BUCKET_SPAN;
 299     // return invalid here instead of clipping, so callers can discard
 300     // invalid buckets without having to check if it's an existing one
 301     if ((clat < -90.0) || (clat > 90.0)) {
 302         return 0;
 303     }
 304
 305     // find the lon span for the new latitude
 306     double trg_span = sg_bucket_span( clat );
 307
 308     // if target span < src_span, return multiple buckets...
 309     if ( trg_span < src_span ) {
 310         // calc center longitude of westernmost sibling
 311         double start_lon = get_center_lat() - src_span/2 + trg_span/2;
 312
 313         unsigned int num_buckets = src_span/trg_span;
 314         for ( unsigned int x = 0; x < num_buckets; x++ ) {
 315             double tmp = start_lon + x * trg_span;
 316             tmp = SGMiscd::normalizePeriodic(-180.0, 180.0, tmp);
 317
 318             SGBucket b;
 319             b.innerSet(tmp, clat);
 320
 321             buckets.push_back( b );
 322         }
 323     } else {
 324         // just return the single sibling
 325         double tmp = get_center_lon() + dx * trg_span;
 326         tmp = SGMiscd::normalizePeriodic(-180.0, 180.0, tmp);
 327
 328         SGBucket b;
 329         b.innerSet(tmp, clat);
 330
 331         buckets.push_back( b );
 332     }
 333
 334     return buckets.size();
 335 }
 336
 337 SGBucket SGBucket::sibling(int dx, int dy) const
 338 {
 339     if (!isValid()) {
 340         SG_LOG(SG_TERRAIN, SG_WARN, "SGBucket::sibling: requesting sibling of invalid bucket");
 341         return SGBucket();
 342     }
 343
 344     double clat = get_center_lat() + dy * SG_BUCKET_SPAN;
 345     // return invalid here instead of clipping, so callers can discard
 346     // invalid buckets without having to check if it's an existing one
 347     if ((clat < -90.0) || (clat > 90.0)) {
 348         return SGBucket();
 349     }
 350
 351     // find the lon span for the new latitude
 352     double span = sg_bucket_span( clat );
 353
 354     double tmp = get_center_lon() + dx * span;
 355     tmp = SGMiscd::normalizePeriodic(-180.0, 180.0, tmp);
 356
 357     SGBucket b;
 358     b.innerSet(tmp, clat);
 359     return b;
 360 }
 361
 362 std::string SGBucket::gen_index_str() const
 363 {
 364         char tmp[20];
 365         ::snprintf(tmp, 20, "%ld",
 366                  (((long)lon + 180) << 14) + ((lat + 90) << 6)
 367                  + (y << 3) + x);
 368         return (std::string)tmp;
 369 }
 370
 371 #ifndef NO_DEPRECATED_API
 372 // find the bucket which is offset by the specified tile units in the
 373 // X & Y direction.  We need the current lon and lat to resolve
 374 // ambiguities when going from a wider tile to a narrower one above or
 375 // below.  This assumes that we are feeding in
 376 SGBucket sgBucketOffset( double dlon, double dlat, int dx, int dy ) {
 377     SGBucket result( dlon, dlat );
 378     double clat = result.get_center_lat() + dy * SG_BUCKET_SPAN;
 379
 380     // walk dy units in the lat direction
 381     result.set_bucket( dlon, clat );
 382
 383     // find the lon span for the new latitude
 384     double span = sg_bucket_span( clat );
 385
 386     // walk dx units in the lon direction
 387     double tmp = dlon + dx * span;
 388     while ( tmp < -180.0 ) {
 389         tmp += 360.0;
 390     }
 391     while ( tmp >= 180.0 ) {
 392         tmp -= 360.0;
 393     }
 394     result.set_bucket( tmp, clat );
 395
 396     return result;
 397 }
 398 #endif
 399
 400 // calculate the offset between two buckets
 401 void sgBucketDiff( const SGBucket& b1, const SGBucket& b2, int *dx, int *dy ) {
 402
 403     // Latitude difference
 404     double c1_lat = b1.get_center_lat();
 405     double c2_lat = b2.get_center_lat();
 406     double diff_lat = c2_lat - c1_lat;
 407
 408 #ifdef HAVE_RINT
 409     *dy = (int)rint( diff_lat / SG_BUCKET_SPAN );
 410 #else
 411     if ( diff_lat > 0 ) {
 412         *dy = (int)( diff_lat / SG_BUCKET_SPAN + 0.5 );
 413     } else {
 414         *dy = (int)( diff_lat / SG_BUCKET_SPAN - 0.5 );
 415     }
 416 #endif
 417
 418     // longitude difference
 419     double diff_lon=0.0;
 420     double span=0.0;
 421
 422     SGBucket tmp_bucket;
 423     // To handle crossing the bucket size boundary
 424     //  we need to account for different size buckets.
 425
 426     if ( sg_bucket_span(c1_lat) <= sg_bucket_span(c2_lat) )
 427     {
 428         span = sg_bucket_span(c1_lat);
 429     } else {
 430         span = sg_bucket_span(c2_lat);
 431     }
 432
 433     diff_lon = b2.get_center_lon() - b1.get_center_lon();
 434
 435     if (diff_lon <0.0)
 436     {
 437        diff_lon -= b1.get_width()*0.5 + b2.get_width()*0.5 - span;
 438     }
 439     else
 440     {
 441        diff_lon += b1.get_width()*0.5 + b2.get_width()*0.5 - span;
 442     }
 443
 444
 445 #ifdef HAVE_RINT
 446     *dx = (int)rint( diff_lon / span );
 447 #else
 448     if ( diff_lon > 0 ) {
 449         *dx = (int)( diff_lon / span + 0.5 );
 450     } else {
 451         *dx = (int)( diff_lon / span - 0.5 );
 452     }
 453 #endif
 454 }
 455
 456 void sgGetBuckets( const SGGeod& min, const SGGeod& max, std::vector<SGBucket>& list ) {
 457     double lon, lat, span;
 458
 459     for (lat = min.getLatitudeDeg(); lat < max.getLatitudeDeg()+SG_BUCKET_SPAN; lat += SG_BUCKET_SPAN) {
 460         span = sg_bucket_span( lat );
 461         for (lon = min.getLongitudeDeg(); lon <= max.getLongitudeDeg(); lon += span)
 462         {
 463             SGBucket b(SGGeod::fromDeg(lon, lat));
 464             if (!b.isValid()) {
 465                 continue;
 466             }
 467
 468             list.push_back(b);
 469         }
 470     }
 471 }
 472
 473 std::ostream& operator<< ( std::ostream& out, const SGBucket& b )
 474 {
 475     return out << b.lon << ":" << (int)b.x << ", " << b.lat << ":" << (int)b.y;
 476 }
 477