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