gov/nasa/worldwind/globes/FlatGlobe.java

   1 /*
   2 Copyright (C) 2001, 2006 United States Government
   3 as represented by the Administrator of the
   4 National Aeronautics and Space Administration.
   5 All Rights Reserved.
   6 */
   7 package gov.nasa.worldwind.globes;
   8
   9 import gov.nasa.worldwind.*;
  10 import gov.nasa.worldwind.avlist.AVKey;
  11 import gov.nasa.worldwind.render.DrawContext;
  12 import gov.nasa.worldwind.geom.*;
  13 import gov.nasa.worldwind.util.Logging;
  14
  15 /** Experimental flat globe
  16  * See TODOs for difference from EllipsoidalGlobe
  17  * @author Patrick Murris - base on EllipsoidalGlobe
  18  * @version $Id$
  19  */
  20 public class FlatGlobe extends WWObjectImpl implements Globe
  21 {
  22     public final static String PROJECTION_LAT_LON = "gov.nasa.worldwind.globes.projectionLatLon";
  23     public final static String PROJECTION_SINUSOIDAL = "gov.nasa.worldwind.globes.projectionSinusoidal";
  24     public final static String PROJECTION_MERCATOR = "gov.nasa.worldwind.globes.projectionMercator";
  25     public final static String PROJECTION_TEST = "gov.nasa.worldwind.globes.projectionTest";
  26
  27     private final double equatorialRadius;
  28     private final double polarRadius;
  29     private final double es;
  30     private final Vec4 center;
  31     private Tessellator tessellator;
  32
  33     private final ElevationModel elevationModel;
  34     private String projection = PROJECTION_LAT_LON;
  35
  36     public FlatGlobe(double equatorialRadius, double polarRadius, double es, ElevationModel em)
  37     {
  38         if (em == null)
  39         {
  40             String msg = Logging.getMessage("nullValue.ElevationModelIsNull");
  41             Logging.logger().severe(msg);
  42             throw new IllegalArgumentException(msg);
  43         }
  44
  45         this.equatorialRadius = equatorialRadius;
  46         this.polarRadius = polarRadius;
  47         this.es = es; // assume it's consistent with the two radii
  48         this.center = Vec4.ZERO;
  49         this.elevationModel = em;
  50     }
  51
  52     public Tessellator getTessellator()
  53     {
  54         return tessellator;
  55     }
  56
  57     public void setTessellator(Tessellator tessellator)
  58     {
  59         this.tessellator = tessellator;
  60     }
  61
  62     public final double getRadius()
  63     {
  64         return this.equatorialRadius;
  65     }
  66
  67     public final double getEquatorialRadius()
  68     {
  69         return this.equatorialRadius;
  70     }
  71
  72     public final double getPolarRadius()
  73     {
  74         return this.polarRadius;
  75     }
  76
  77     public double getMaximumRadius()
  78     {
  79         return this.equatorialRadius;
  80     }
  81
  82     // TODO: Find a more accurate workaround then getRadius()
  83     public double getRadiusAt(Angle latitude, Angle longitude)
  84     {
  85         if (latitude == null || longitude == null)
  86         {
  87             String msg = Logging.getMessage("nullValue.AngleIsNull");
  88             Logging.logger().severe(msg);
  89             throw new IllegalArgumentException(msg);
  90         }
  91         return getRadius();
  92         //return this.computePointFromPosition(latitude, longitude, 0d).getLength3();
  93     }
  94
  95     // TODO: Find a more accurate workaround then getRadius()
  96     public double getRadiusAt(LatLon latLon)
  97     {
  98         if (latLon == null)
  99         {
 100             String msg = Logging.getMessage("nullValue.LatLonIsNull");
 101             Logging.logger().severe(msg);
 102             throw new IllegalArgumentException(msg);
 103         }
 104         return getRadius();
 105         //return this.computePointFromPosition(latLon.getLatitude(), latLon.getLongitude(), 0d).getLength3();
 106     }
 107
 108     public double getEccentricitySquared()
 109     {
 110         return this.es;
 111     }
 112
 113     public final double getDiameter()
 114     {
 115         return this.equatorialRadius * 2;
 116     }
 117
 118     public final Vec4 getCenter()
 119     {
 120         return this.center;
 121     }
 122
 123     public double getMaxElevation()
 124     {
 125         return this.elevationModel.getMaxElevation();
 126     }
 127
 128     public double getMinElevation()
 129     {
 130         return this.elevationModel.getMinElevation();
 131     }
 132
 133     public double getMaxElevation(Sector sector)
 134     {
 135         return this.elevationModel.getMaxElevation(sector);
 136     }
 137
 138     public double getMinElevation(Sector sector)
 139     {
 140         return this.elevationModel.getMinElevation(sector);
 141     }
 142
 143     public final Extent getExtent()
 144     {
 145         return this;
 146     }
 147
 148     public void setProjection(String projection)
 149     {
 150         this.projection = projection;
 151     }
 152
 153     public String getProjection()
 154     {
 155         return this.projection;
 156     }
 157
 158     public boolean intersects(Frustum frustum)
 159     {
 160         if (frustum == null)
 161         {
 162             String message = Logging.getMessage("nullValue.FrustumIsNull");
 163             Logging.logger().severe(message);
 164             throw new IllegalArgumentException(message);
 165         }
 166
 167         return frustum.intersects(this);
 168     }
 169
 170     public Intersection[] intersect(Line line)
 171     {
 172         return this.intersect(line, this.equatorialRadius, this.polarRadius);
 173     }
 174
 175     public Intersection[] intersect(Line line, double altitude)
 176     {
 177         return this.intersect(line, this.equatorialRadius + altitude, this.polarRadius + altitude);
 178     }
 179
 180     // TODO: plane/line intersection point (OK)
 181     // TODO: extract altitude from equRadius by subtracting this.equatorialRadius (OK)
 182     private Intersection[] intersect(Line line, double equRadius, double polRadius)
 183     {
 184         if (line == null)
 185         {
 186             String message = Logging.getMessage("nullValue.LineIsNull");
 187             Logging.logger().severe(message);
 188             throw new IllegalArgumentException(message);
 189         }
 190         // Intersection with world plane
 191         Plane plane = new Plane(0, 0, 1, equRadius - this.equatorialRadius);   // Flat globe plane
 192         Vec4 p = plane.intersect(line);
 193         if(p == null)
 194             return null;
 195         // Check if we are in the world boundaries
 196         Position pos = this.computePositionFromPoint(p);
 197         if (pos == null)
 198             return null;
 199         if(pos.getLatitude().degrees < -90 || pos.getLatitude().degrees > 90 ||
 200                 pos.getLongitude().degrees < -180 || pos.getLongitude().degrees > 180)
 201             return null;
 202
 203         return new Intersection[] {new Intersection(p, false)};
 204     }
 205
 206     static private double discriminant(double a, double b, double c)
 207     {
 208         return b * b - 4 * a * c;
 209     }
 210
 211     // TODO: plane/line intersection test (OK)
 212     public boolean intersects(Line line)
 213     {
 214         if (line == null)
 215         {
 216             String msg = Logging.getMessage("nullValue.LineIsNull");
 217             Logging.logger().severe(msg);
 218             throw new IllegalArgumentException(msg);
 219         }
 220
 221         return this.intersect(line) != null;
 222     }
 223
 224     // TODO: plane/plane intersection test (OK)
 225     public boolean intersects(Plane plane)
 226     {
 227         if (plane == null)
 228         {
 229             String msg = Logging.getMessage("nullValue.PlaneIsNull");
 230             Logging.logger().severe(msg);
 231             throw new IllegalArgumentException(msg);
 232         }
 233
 234         return !plane.getNormal().equals(Vec4.UNIT_Z);
 235
 236     }
 237
 238     // TODO: return constant (OK)
 239     public Vec4 computeSurfaceNormalAtPoint(Vec4 p)
 240     {
 241         if (p == null)
 242         {
 243             String msg = Logging.getMessage("nullValue.PointIsNull");
 244             Logging.logger().severe(msg);
 245             throw new IllegalArgumentException(msg);
 246         }
 247
 248         return new Vec4(0, 0, 1);
 249     }
 250
 251     public final ElevationModel getElevationModel()
 252     {
 253         return this.elevationModel;
 254     }
 255
 256     public Double getBestElevation(Angle latitude, Angle longitude)
 257     {
 258         return this.elevationModel.getBestElevation(latitude, longitude);
 259     }
 260
 261     // TODO: return zero if outside the lat/lon normal boundaries (OK)
 262     public final Double getElevationAtResolution(Angle latitude, Angle longitude, double resolution)
 263     {
 264         if (latitude == null || longitude == null)
 265         {
 266             String message = Logging.getMessage("nullValue.LatitudeOrLongitudeIsNull");
 267             Logging.logger().severe(message);
 268             throw new IllegalArgumentException(message);
 269         }
 270         if(latitude.degrees < -90 || latitude.degrees > 90 || longitude.degrees < -180 || longitude.degrees > 180)
 271             return null;
 272
 273         int target = this.elevationModel.getTargetResolution(this, resolution);
 274         return this.elevationModel.getElevationAtResolution(latitude, longitude, target);
 275     }
 276
 277     // TODO: return zero if outside the lat/lon normal boundaries (OK)
 278     public final double getElevation(Angle latitude, Angle longitude)
 279     {
 280         if (latitude == null || longitude == null)
 281         {
 282             String message = Logging.getMessage("nullValue.LatitudeOrLongitudeIsNull");
 283             Logging.logger().severe(message);
 284             throw new IllegalArgumentException(message);
 285         }
 286         if(latitude.degrees < -90 || latitude.degrees > 90 || longitude.degrees < -180 || longitude.degrees > 180)
 287             return 0;
 288         return this.elevationModel != null ? this.elevationModel.getElevation(latitude, longitude) : 0;
 289     }
 290
 291     public final Vec4 computePointFromPosition(Position position)
 292     {
 293         if (position == null)
 294         {
 295             String message = Logging.getMessage("nullValue.PositionIsNull");
 296             Logging.logger().severe(message);
 297             throw new IllegalArgumentException(message);
 298         }
 299
 300         return this.geodeticToCartesian(position.getLatitude(), position.getLongitude(), position.getElevation());
 301     }
 302
 303     public final Vec4 computePointFromPosition(Angle latitude, Angle longitude, double metersElevation)
 304     {
 305         if (latitude == null || longitude == null)
 306         {
 307             String message = Logging.getMessage("nullValue.LatitudeOrLongitudeIsNull");
 308             Logging.logger().severe(message);
 309             throw new IllegalArgumentException(message);
 310         }
 311
 312         return this.geodeticToCartesian(latitude, longitude, metersElevation);
 313     }
 314
 315     public final Position computePositionFromPoint(Vec4 point)
 316     {
 317         if (point == null)
 318         {
 319             String message = Logging.getMessage("nullValue.PointIsNull");
 320             Logging.logger().severe(message);
 321             throw new IllegalArgumentException(message);
 322         }
 323
 324         return this.cartesianToGeodetic(point);
 325     }
 326
 327     public final Position getIntersectionPosition(Line line)
 328     {
 329         if (line == null)
 330         {
 331             String msg = Logging.getMessage("nullValue.LineIsNull");
 332             Logging.logger().severe(msg);
 333             throw new IllegalArgumentException(msg);
 334         }
 335
 336         Intersection[] intersections = this.intersect(line);
 337         if (intersections == null)
 338             return null;
 339
 340         return this.computePositionFromPoint(intersections[0].getIntersectionPoint());
 341     }
 342
 343     /*
 344     // The code below maps latitude / longitude position to globe-centered Cartesian coordinates.
 345     // The Y axis points to the north pole. The Z axis points to the intersection of the prime
 346     // meridian and the equator, in the equatorial plane. The X axis completes a right-handed
 347     // coordinate system, and is 90 degrees east of the Z axis and also in the equatorial plane.
 348
 349     private Vec4 geodeticToCartesianEl(Angle latitude, Angle longitude, double metersElevation)
 350     {
 351         if (latitude == null || longitude == null)
 352         {
 353             String message = Logging.getMessage("nullValue.LatitudeOrLongitudeIsNull");
 354             Logging.logger().severe(message);
 355             throw new IllegalArgumentException(message);
 356         }
 357
 358         double cosLat = latitude.cos();
 359         double sinLat = latitude.sin();
 360
 361         double rpm = // getRadius (in meters) of vertical in prime meridian
 362             this.equatorialRadius / Math.sqrt(1.0 - this.es * sinLat * sinLat);
 363
 364         double x = (rpm + metersElevation) * cosLat * longitude.sin();
 365         double y = (rpm * (1.0 - this.es) + metersElevation) * sinLat;
 366         double z = (rpm + metersElevation) * cosLat * longitude.cos();
 367
 368         Vec4 cart = new Vec4(x, y, z);
 369         //System.out.println("geodeticToCartesian: " + latitude + ", " + longitude + ", " + metersElevation + " / " + cart);
 370         return cart;
 371     }
 372
 373     private Position cartesianToGeodeticEl(Vec4 cart)
 374     {
 375         if (cart == null)
 376         {
 377             String message = Logging.getMessage("nullValue.PointIsNull");
 378             Logging.logger().severe(message);
 379             throw new IllegalArgumentException(message);
 380         }
 381
 382         // according to
 383         // H. Vermeille,
 384         // Direct transformation from geocentric to geodetic ccordinates,
 385         // Journal of Geodesy (2002) 76:451-454
 386         double ra2 = 1 / (this.equatorialRadius * equatorialRadius);
 387
 388         double X = cart.z;
 389         //noinspection SuspiciousNameCombination
 390         double Y = cart.x;
 391         double Z = cart.y;
 392         double e2 = this.es;
 393         double e4 = e2 * e2;
 394
 395         double XXpYY = X * X + Y * Y;
 396         double sqrtXXpYY = Math.sqrt(XXpYY);
 397         double p = XXpYY * ra2;
 398         double q = Z * Z * (1 - e2) * ra2;
 399         double r = 1 / 6.0 * (p + q - e4);
 400         double s = e4 * p * q / (4 * r * r * r);
 401         double t = Math.pow(1 + s + Math.sqrt(s * (2 + s)), 1 / 3.0);
 402         double u = r * (1 + t + 1 / t);
 403         double v = Math.sqrt(u * u + e4 * q);
 404         double w = e2 * (u + v - q) / (2 * v);
 405         double k = Math.sqrt(u + v + w * w) - w;
 406         double D = k * sqrtXXpYY / (k + e2);
 407         double lon = 2 * Math.atan2(Y, X + sqrtXXpYY);
 408         double sqrtDDpZZ = Math.sqrt(D * D + Z * Z);
 409         double lat = 2 * Math.atan2(Z, D + sqrtDDpZZ);
 410         double elevation = (k + e2 - 1) * sqrtDDpZZ / k;
 411
 412         Position pos = Position.fromRadians(lat, lon, elevation);
 413         //System.out.println("cartesianToGeodetic: " + cart + " / " + pos);
 414         return pos;
 415     }  */
 416
 417     // The code below maps latitude / longitude position to globe-centered Cartesian coordinates.
 418     // The Y axis points to the north pole. The Z axis points to the intersection of the prime
 419     // meridian and the equator, in the equatorial plane. The X axis completes a right-handed
 420     // coordinate system, and is 90 degrees east of the Z axis and also in the equatorial plane.
 421     // TODO: Implement flat projections (OK)
 422     private Vec4 geodeticToCartesian(Angle latitude, Angle longitude, double metersElevation)
 423     {
 424         if (latitude == null || longitude == null)
 425         {
 426             String message = Logging.getMessage("nullValue.LatitudeOrLongitudeIsNull");
 427             Logging.logger().severe(message);
 428             throw new IllegalArgumentException(message);
 429         }
 430
 431         Vec4 cart = null;
 432         if(this.projection.compareToIgnoreCase(PROJECTION_LAT_LON) == 0)
 433         {
 434             // Lat/Lon projection - plate carré
 435             cart = new Vec4(this.equatorialRadius * longitude.radians,
 436                     this.equatorialRadius * latitude.radians,
 437                     metersElevation);
 438         }
 439         else if(this.projection.compareToIgnoreCase(PROJECTION_MERCATOR) == 0)
 440         {
 441             // Mercator projection
 442             if(latitude.degrees > 75) latitude = Angle.fromDegrees(75);
 443             if(latitude.degrees < -75) latitude = Angle.fromDegrees(-75);
 444             cart = new Vec4(this.equatorialRadius * longitude.radians,
 445                     this.equatorialRadius * Math.log(Math.tan(Math.PI / 4 + latitude.radians / 2)),
 446                     metersElevation);
 447         }
 448         else if(this.projection.compareToIgnoreCase(PROJECTION_SINUSOIDAL) == 0)
 449         {
 450             // Sinusoidal projection
 451             cart = new Vec4(this.equatorialRadius * longitude.radians * latitude.cos(),
 452                     this.equatorialRadius * latitude.radians,
 453                     metersElevation);
 454         }
 455         else if(this.projection.compareToIgnoreCase(PROJECTION_TEST) == 0)
 456         {
 457             // Test projection
 458             double r = Math.sqrt(latitude.radians * latitude.radians + longitude.radians * longitude.radians) * this.equatorialRadius;
 459             cart = new Vec4(this.equatorialRadius * longitude.radians * Math.pow(latitude.cos(), .3),
 460                     this.equatorialRadius * latitude.radians,
 461                     metersElevation);
 462         }
 463         //System.out.println("geodeticToCartesian: " + latitude + ", " + longitude + ", " + metersElevation + " / " + cart);
 464         //System.out.println(this.projection);
 465         return cart;
 466     }
 467
 468     // TODO: Implement flat projections (OK)
 469     private Position cartesianToGeodetic(Vec4 cart)
 470     {
 471         if (cart == null)
 472         {
 473             String message = Logging.getMessage("nullValue.PointIsNull");
 474             Logging.logger().severe(message);
 475             throw new IllegalArgumentException(message);
 476         }
 477
 478         Position pos = null;
 479         if(this.projection.compareToIgnoreCase(PROJECTION_LAT_LON) == 0)
 480         {
 481         // Lat/Lon projection - plate carré
 482         pos =  Position.fromRadians(
 483                 cart.y / this.equatorialRadius,
 484                 cart.x / this.equatorialRadius,
 485                 cart.z);
 486         }
 487         else if(this.projection.compareToIgnoreCase(PROJECTION_MERCATOR) == 0)
 488         {
 489         // Mercator projection
 490         pos =  Position.fromRadians(
 491                 Math.atan(Math.sinh(cart.y / this.equatorialRadius)) ,
 492                 cart.x / this.equatorialRadius,
 493                 cart.z);
 494         }
 495         else if(this.projection.compareToIgnoreCase(PROJECTION_SINUSOIDAL) == 0)
 496         {
 497         // Sinusoidal projection
 498         pos =  Position.fromRadians(
 499                 cart.y / this.equatorialRadius,
 500                 cart.x / this.equatorialRadius / Angle.fromRadians(cart.y / this.equatorialRadius).cos(),
 501                 cart.z);
 502         }
 503         else if(this.projection.compareToIgnoreCase(PROJECTION_TEST) == 0)
 504         {
 505         // Test projection
 506         pos =  Position.fromRadians(
 507                 cart.y / this.equatorialRadius,
 508                 cart.x / this.equatorialRadius / Math.pow(Angle.fromRadians(cart.y / this.equatorialRadius).cos(), .3),
 509                 cart.z);
 510         }
 511         //System.out.println("cartesianToGeodetic: " + cart + " / " + pos);
 512         return pos;
 513     }
 514
 515
 516     /**
 517      * Returns a cylinder that minimally surrounds the sector at a specified vertical exaggeration.
 518      *
 519      * @param verticalExaggeration the vertical exaggeration to apply to the globe's elevations when computing the
 520      *                             cylinder.
 521      * @param sector               the sector to return the bounding cylinder for.
 522      * @return The minimal bounding cylinder in Cartesian coordinates.
 523      * @throws IllegalArgumentException if <code>globe</code> or <code>sector</code> is null
 524      */
 525     // TODO: Adapt to flat world tiles (OK)
 526     public Cylinder computeBoundingCylinder(double verticalExaggeration, Sector sector)
 527     {
 528         if (sector == null)
 529         {
 530             String msg = Logging.getMessage("nullValue.SectorIsNull");
 531             Logging.logger().severe(msg);
 532             throw new IllegalArgumentException(msg);
 533         }
 534
 535         // Compute the center points of the bounding cylinder's top and bottom planes.
 536         LatLon center = sector.getCentroid();
 537         double maxHeight = this.getMaxElevation(sector) * verticalExaggeration;
 538         double minHeight = 0; //globe.getMinElevation(sector) * verticalExaggeration;
 539         Vec4 centroidTop = this.computePointFromPosition(center.getLatitude(), center.getLongitude(), maxHeight);
 540         Vec4 centroidBot = this.computePointFromPosition(center.getLatitude(), center.getLongitude(), minHeight);
 541
 542         // Compute radius of circumscribing circle around general quadrilateral.
 543         Vec4 northwest = this.computePointFromPosition(sector.getMaxLatitude(), sector.getMinLongitude(), maxHeight);
 544         Vec4 southeast = this.computePointFromPosition(sector.getMinLatitude(), sector.getMaxLongitude(), maxHeight);
 545         Vec4 southwest = this.computePointFromPosition(sector.getMinLatitude(), sector.getMinLongitude(), maxHeight);
 546         Vec4 northeast = this.computePointFromPosition(sector.getMaxLatitude(), sector.getMaxLongitude(), maxHeight);
 547         double a = southwest.distanceTo3(southeast);
 548         double b = southeast.distanceTo3(northeast);
 549         double c = northeast.distanceTo3(northwest);
 550         double d = northwest.distanceTo3(southwest);
 551         double s = 0.5 * (a + b + c + d);
 552         double area = Math.sqrt((s - a) * (s - b) * (s - c) * (s - d));
 553         double radius = Math.sqrt((a * b + c * d) * (a * d + b * c) * (a * c + b * d)) / (4d * area);
 554
 555         return new Cylinder(centroidBot, centroidTop, radius);
 556     }
 557
 558     public boolean equals(Object o)
 559     {
 560         if (this == o) return true;
 561         if (o == null || getClass() != o.getClass()) return false;
 562
 563         FlatGlobe flatGlobe = (FlatGlobe) o;
 564
 565         if (Double.compare(flatGlobe.equatorialRadius, equatorialRadius) != 0) return false;
 566         if (Double.compare(flatGlobe.polarRadius, polarRadius) != 0) return false;
 567         if (center != null ? !center.equals(flatGlobe.center) : flatGlobe.center != null) return false;
 568         if (elevationModel != null ? !elevationModel.equals(flatGlobe.elevationModel)
 569             : flatGlobe.elevationModel != null)
 570             return false;
 571         if (projection != null ? !projection.equals(flatGlobe.projection) : flatGlobe.projection != null) return false;
 572         if (tessellator != null ? !tessellator.equals(flatGlobe.tessellator) : flatGlobe.tessellator != null)
 573             return false;
 574
 575         return true;
 576     }
 577
 578     public int hashCode()
 579     {
 580         int result;
 581         long temp;
 582         temp = equatorialRadius != +0.0d ? Double.doubleToLongBits(equatorialRadius) : 0L;
 583         result = (int) (temp ^ (temp >>> 32));
 584         temp = polarRadius != +0.0d ? Double.doubleToLongBits(polarRadius) : 0L;
 585         result = 31 * result + (int) (temp ^ (temp >>> 32));
 586         result = 31 * result + (center != null ? center.hashCode() : 0);
 587         result = 31 * result + (tessellator != null ? tessellator.hashCode() : 0);
 588         result = 31 * result + (elevationModel != null ? elevationModel.hashCode() : 0);
 589         result = 31 * result + (projection != null ? projection.hashCode() : 0);
 590         return result;
 591     }
 592
 593     public SectorGeometryList tessellate(DrawContext dc)
 594     {
 595         if (this.tessellator == null)
 596         {
 597             this.tessellator = (Tessellator) WorldWind.createConfigurationComponent(AVKey.TESSELLATOR_CLASS_NAME);
 598         }
 599
 600         return this.tessellator.tessellate(dc);
 601     }
 602 }