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 <math.h>
  31
  32 #include <simgear/misc/sg_path.hxx>
  33
  34 #include "newbucket.hxx"
  35
  36
  37 // default constructor
  38 SGBucket::SGBucket() {
  39 }
  40
  41
  42 // constructor for specified location
  43 SGBucket::SGBucket(const double dlon, const double dlat) {
  44     set_bucket(dlon, dlat);
  45 }
  46
  47 SGBucket::SGBucket(const SGGeod& geod) {
  48     set_bucket(geod);
  49 }
  50
  51 // create an impossible bucket if false
  52 SGBucket::SGBucket(const bool is_good) {
  53     set_bucket(0.0, 0.0);
  54     if ( !is_good ) {
  55         lon = -1000;
  56     }
  57 }
  58
  59
  60 // Parse a unique scenery tile index and find the lon, lat, x, and y
  61 SGBucket::SGBucket(const long int bindex) {
  62     long int index = bindex;
  63
  64     lon = index >> 14;
  65     index -= lon << 14;
  66     lon -= 180;
  67
  68     lat = index >> 6;
  69     index -= lat << 6;
  70     lat -= 90;
  71
  72     y = index >> 3;
  73     index -= y << 3;
  74
  75     x = index;
  76 }
  77
  78
  79 // Set the bucket params for the specified lat and lon
  80 void SGBucket::set_bucket( double dlon, double dlat ) {
  81     //
  82     // latitude first
  83     //
  84     double span = sg_bucket_span( dlat );
  85     double diff = dlon - (double)(int)dlon;
  86
  87     // cout << "diff = " << diff << "  span = " << span << endl;
  88
  89     /* Calculate the greatest integral longitude less than
  90      * or equal to the given longitude (floor(dlon)),
  91      * but attribute coordinates near the east border
  92      * to the next tile.
  93      */
  94     if ( (dlon >= 0) || (fabs(diff) < SG_EPSILON) ) {
  95         lon = (int)dlon;
  96     } else {
  97         lon = (int)dlon - 1;
  98     }
  99
 100     // find subdivision or super lon if needed
 101     if ( span < SG_EPSILON ) {
 102         /* sg_bucket_span() never returns 0.0
 103          * or anything near it, so this really
 104          * should not occur at any time.
 105          */
 106         // polar cap
 107         lon = 0;
 108         x = 0;
 109     } else if ( span <= 1.0 ) {
 110         /* We have more than one tile per degree of
 111          * longitude, so we need an x offset.
 112          */
 113         x = (int)((dlon - lon) / span);
 114     } else {
 115         /* We have one or more degrees per tile,
 116          * so we need to find the base longitude
 117          * of that tile.
 118          *
 119          * First we calculate the integral base longitude
 120          * (e.g. -85.5 => -86) and then find the greatest
 121          * multiple of span that is less than or equal to
 122          * that longitude.
 123          *
 124          * That way, the Greenwich Meridian is always
 125          * a tile border.
 126          *
 127          * This gets us into trouble with the polar caps,
 128          * which have width 360 and thus either span
 129          * the range from 0 to 360 or from -360 to 0
 130          * degrees, depending on whether lon is positive
 131          * or negative!
 132          *
 133          * We also get into trouble with the 8 degree tiles
 134          * north of 88N and south of 88S, because the west-
 135          * and east-most tiles in that range will cover 184W
 136          * to 176W and 176E to 184E respectively, with their
 137          * center at 180E/W!
 138          */
 139         lon=(int)floor(floor((lon+SG_EPSILON)/span)*span);
 140         /* Correct the polar cap issue */
 141         if ( lon < -180 ) {
 142             lon = -180;
 143         }
 144         x = 0;
 145     }
 146
 147     //
 148     // then latitude
 149     //
 150     diff = dlat - (double)(int)dlat;
 151
 152     /* Again, a modified floor() function (see longitude) */
 153     if ( (dlat >= 0) || (fabs(diff) < SG_EPSILON) ) {
 154         lat = (int)dlat;
 155     } else {
 156         lat = (int)dlat - 1;
 157     }
 158     /* Latitude base and offset are easier, as
 159      * tiles always are 1/8 degree of latitude wide.
 160      */
 161     y = (int)((dlat - lat) * 8);
 162 }
 163
 164
 165 void SGBucket::set_bucket(const SGGeod& geod)
 166 {
 167     set_bucket(geod.getLongitudeDeg(), geod.getLatitudeDeg());
 168 }
 169
 170 // Build the path name for this bucket
 171 std::string SGBucket::gen_base_path() const {
 172     // long int index;
 173     int top_lon, top_lat, main_lon, main_lat;
 174     char hem, pole;
 175     char raw_path[256];
 176
 177     top_lon = lon / 10;
 178     main_lon = lon;
 179     if ( (lon < 0) && (top_lon * 10 != lon) ) {
 180         top_lon -= 1;
 181     }
 182     top_lon *= 10;
 183     if ( top_lon >= 0 ) {
 184         hem = 'e';
 185     } else {
 186         hem = 'w';
 187         top_lon *= -1;
 188     }
 189     if ( main_lon < 0 ) {
 190         main_lon *= -1;
 191     }
 192
 193     top_lat = lat / 10;
 194     main_lat = lat;
 195     if ( (lat < 0) && (top_lat * 10 != lat) ) {
 196         top_lat -= 1;
 197     }
 198     top_lat *= 10;
 199     if ( top_lat >= 0 ) {
 200         pole = 'n';
 201     } else {
 202         pole = 's';
 203         top_lat *= -1;
 204     }
 205     if ( main_lat < 0 ) {
 206         main_lat *= -1;
 207     }
 208
 209     snprintf(raw_path, 256, "%c%03d%c%02d/%c%03d%c%02d",
 210             hem, top_lon, pole, top_lat,
 211             hem, main_lon, pole, main_lat);
 212
 213     SGPath path( raw_path );
 214
 215     return path.str();
 216 }
 217
 218
 219 // return width of the tile in degrees
 220 double SGBucket::get_width() const {
 221     return sg_bucket_span( get_center_lat() );
 222 }
 223
 224
 225 // return height of the tile in degrees
 226 double SGBucket::get_height() const {
 227     return SG_BUCKET_SPAN;
 228 }
 229
 230
 231 // return width of the tile in meters
 232 double SGBucket::get_width_m() const {
 233     double clat = (int)get_center_lat();
 234     if ( clat > 0 ) {
 235         clat = (int)clat + 0.5;
 236     } else {
 237         clat = (int)clat - 0.5;
 238     }
 239     double clat_rad = clat * SGD_DEGREES_TO_RADIANS;
 240     double cos_lat = cos( clat_rad );
 241     double local_radius = cos_lat * SG_EQUATORIAL_RADIUS_M;
 242     double local_perimeter = local_radius * SGD_2PI;
 243     double degree_width = local_perimeter / 360.0;
 244
 245     return get_width() * degree_width;
 246 }
 247
 248
 249 // return height of the tile in meters
 250 double SGBucket::get_height_m() const {
 251     double perimeter = SG_EQUATORIAL_RADIUS_M * SGD_2PI;
 252     double degree_height = perimeter / 360.0;
 253
 254     return SG_BUCKET_SPAN * degree_height;
 255 }
 256
 257
 258 // find the bucket which is offset by the specified tile units in the
 259 // X & Y direction.  We need the current lon and lat to resolve
 260 // ambiguities when going from a wider tile to a narrower one above or
 261 // below.  This assumes that we are feeding in
 262 SGBucket sgBucketOffset( double dlon, double dlat, int dx, int dy ) {
 263     SGBucket result( dlon, dlat );
 264     double clat = result.get_center_lat() + dy * SG_BUCKET_SPAN;
 265
 266     // walk dy units in the lat direction
 267     result.set_bucket( dlon, clat );
 268
 269     // find the lon span for the new latitude
 270     double span = sg_bucket_span( clat );
 271
 272     // walk dx units in the lon direction
 273     double tmp = dlon + dx * span;
 274     while ( tmp < -180.0 ) {
 275         tmp += 360.0;
 276     }
 277     while ( tmp >= 180.0 ) {
 278         tmp -= 360.0;
 279     }
 280     result.set_bucket( tmp, clat );
 281
 282     return result;
 283 }
 284
 285
 286 // calculate the offset between two buckets
 287 void sgBucketDiff( const SGBucket& b1, const SGBucket& b2, int *dx, int *dy ) {
 288
 289     // Latitude difference
 290     double c1_lat = b1.get_center_lat();
 291     double c2_lat = b2.get_center_lat();
 292     double diff_lat = c2_lat - c1_lat;
 293
 294 #ifdef HAVE_RINT
 295     *dy = (int)rint( diff_lat / SG_BUCKET_SPAN );
 296 #else
 297     if ( diff_lat > 0 ) {
 298         *dy = (int)( diff_lat / SG_BUCKET_SPAN + 0.5 );
 299     } else {
 300         *dy = (int)( diff_lat / SG_BUCKET_SPAN - 0.5 );
 301     }
 302 #endif
 303
 304     // longitude difference
 305     double diff_lon=0.0;
 306     double span=0.0;
 307
 308     SGBucket tmp_bucket;
 309     // To handle crossing the bucket size boundary
 310     //  we need to account for different size buckets.
 311
 312     if ( sg_bucket_span(c1_lat) <= sg_bucket_span(c2_lat) )
 313     {
 314         span = sg_bucket_span(c1_lat);
 315     } else {
 316         span = sg_bucket_span(c2_lat);
 317     }
 318
 319     diff_lon = b2.get_center_lon() - b1.get_center_lon();
 320
 321     if (diff_lon <0.0)
 322     {
 323        diff_lon -= b1.get_width()*0.5 + b2.get_width()*0.5 - span;
 324     }
 325     else
 326     {
 327        diff_lon += b1.get_width()*0.5 + b2.get_width()*0.5 - span;
 328     }
 329
 330
 331 #ifdef HAVE_RINT
 332     *dx = (int)rint( diff_lon / span );
 333 #else
 334     if ( diff_lon > 0 ) {
 335         *dx = (int)( diff_lon / span + 0.5 );
 336     } else {
 337         *dx = (int)( diff_lon / span - 0.5 );
 338     }
 339 #endif
 340 }
 341
 342 void sgGetBuckets( const SGGeod& min, const SGGeod& max, std::vector<SGBucket>& list ) {
 343     double lon, lat, span;
 344
 345     for (lat = min.getLatitudeDeg(); lat <= max.getLatitudeDeg(); lat += SG_BUCKET_SPAN) {
 346         span = sg_bucket_span( lat );
 347         for (lon = min.getLongitudeDeg(); lon <= max.getLongitudeDeg(); lon += span) {
 348             list.push_back( SGBucket(lon , lat) );
 349         }
 350     }
 351 }