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