diff --git a/CHANGES.md b/CHANGES.md index 9b13a5368..9369bc6ac 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -12,6 +12,7 @@ - Added conversions from `std::string` to other metadata types in `MetadataConversions`. This enables the same conversions as `std::string_view`, while allowing runtime engines to use `std::string` for convenience. - Added `applyTextureTransform` property to `TilesetOptions`, which indicates whether to preemptively apply transforms to texture coordinates for textures with the `KHR_texture_transform` extension. - Added `loadGltf` method to `GltfReader`, making it easier to do a full, asynchronous load of a glTF. +- Added `GlobeFlightPath` class to help with calculating fly-to paths. ##### Fixes :wrench: diff --git a/CesiumGeospatial/include/CesiumGeospatial/Ellipsoid.h b/CesiumGeospatial/include/CesiumGeospatial/Ellipsoid.h index 2c4764f35..6d038f2a5 100644 --- a/CesiumGeospatial/include/CesiumGeospatial/Ellipsoid.h +++ b/CesiumGeospatial/include/CesiumGeospatial/Ellipsoid.h @@ -119,6 +119,17 @@ class CESIUMGEOSPATIAL_API Ellipsoid final { std::optional scaleToGeodeticSurface(const glm::dvec3& cartesian) const noexcept; + /** + * @brief Scales the provided cartesian position along the geocentric + * surface normal so that it is on the surface of this ellipsoid. + * + * @param cartesian The cartesian position to scale. + * @retun The scaled position, or the empty optional if the cartesian is at + * the center of this ellipsoid. + */ + std::optional + scaleToGeocentricSurface(const glm::dvec3& cartesian) const noexcept; + /** * @brief The maximum radius in any dimension. * diff --git a/CesiumGeospatial/include/CesiumGeospatial/SimplePlanarEllipsoidCurve.h b/CesiumGeospatial/include/CesiumGeospatial/SimplePlanarEllipsoidCurve.h new file mode 100644 index 000000000..8f9853fb2 --- /dev/null +++ b/CesiumGeospatial/include/CesiumGeospatial/SimplePlanarEllipsoidCurve.h @@ -0,0 +1,102 @@ +#pragma once + +#include "Library.h" + +#include + +#include + +#include + +namespace CesiumGeospatial { + +/** + * @brief Produces points on an ellipse that lies on a plane that intersects the + * center of the earth and each of the input coordinates. The height above the + * surface at each point along the curve will be a linear interpolation between + * the source and destination heights. + */ +class CESIUMGEOSPATIAL_API SimplePlanarEllipsoidCurve final { +public: + /** + * @brief Creates a new instance of {@link SimplePlanarEllipsoidCurve} from a + * source and destination specified in Earth-Centered, Earth-Fixed + * coordinates. + * + * @param ellipsoid The ellipsoid that the source and destination positions + * are relative to. + * @param sourceEcef The position that the path will begin at in ECEF + * coordinates. + * @param destinationEcef The position that the path will end at in ECEF + * coordinates. + * + * @returns An optional type containing a {@link SimplePlanarEllipsoidCurve} + * object representing the generated path, if possible. If it wasn't possible + * to scale the input coordinates to geodetic surface coordinates on a WGS84 + * ellipsoid, this will return {@link std::nullopt} instead. + */ + static std::optional + fromEarthCenteredEarthFixedCoordinates( + const Ellipsoid& ellipsoid, + const glm::dvec3& sourceEcef, + const glm::dvec3& destinationEcef); + + /** + * @brief Creates a new instance of {@link SimplePlanarEllipsoidCurve} from a + * source and destination specified in cartographic coordinates (Longitude, + * Latitude, and Height). + * + * @param ellipsoid The ellipsoid that these cartographic coordinates are + * from. + * @param sourceLlh The position that the path will begin at in Longitude, + * Latitude, and Height. + * @param destinationLlh The position that the path will end at in Longitude, + * Latitude, and Height. + * + * @returns An optional type containing a {@link SimplePlanarEllipsoidCurve} + * object representing the generated path, if possible. If it wasn't possible + * to scale the input coordinates to geodetic surface coordinates on a WGS84 + * ellipsoid, this will return {@link std::nullopt} instead. + */ + static std::optional fromLongitudeLatitudeHeight( + const Ellipsoid& ellipsoid, + const Cartographic& source, + const Cartographic& destination); + + /** + * @brief Samples the curve at the given percentage of its length. + * + * @param percentage The percentage of the curve's length to sample at, + * where 0 is the beginning and 1 is the end. This value will be clamped to + * the range [0..1]. + * @param additionalHeight The height above the earth at this position will be + * calculated by interpolating between the height at the beginning and end of + * the curve based on the value of \p percentage. This parameter specifies an + * additional offset to add to the height. + * + * @returns The position of the given point on this curve in Earth-Centered + * Earth-Fixed coordinates. + */ + glm::dvec3 + getPosition(double percentage, double additionalHeight = 0.0) const; + +private: + SimplePlanarEllipsoidCurve( + const Ellipsoid& ellipsoid, + const glm::dvec3& scaledSourceEcef, + const glm::dvec3& scaledDestinationEcef, + const glm::dvec3& originalSourceEcef, + const glm::dvec3& originalDestinationEcef); + + double _totalAngle; + double _sourceHeight; + double _destinationHeight; + + Ellipsoid _ellipsoid; + glm::dvec3 _sourceDirection; + glm::dvec3 _rotationAxis; + glm::dvec3 _sourceEcef; + glm::dvec3 _destinationEcef; +}; + +} // namespace CesiumGeospatial diff --git a/CesiumGeospatial/src/Ellipsoid.cpp b/CesiumGeospatial/src/Ellipsoid.cpp index 3ae9d0dac..23b233a97 100644 --- a/CesiumGeospatial/src/Ellipsoid.cpp +++ b/CesiumGeospatial/src/Ellipsoid.cpp @@ -145,4 +145,29 @@ Ellipsoid::scaleToGeodeticSurface(const glm::dvec3& cartesian) const noexcept { positionZ * zMultiplier); } +std::optional Ellipsoid::scaleToGeocentricSurface( + const glm::dvec3& cartesian) const noexcept { + + // If the input cartesian is (0, 0, 0), beta will compute to 1.0 / sqrt(0) = + // +Infinity. Let's consider this a failure. + if (Math::equalsEpsilon(glm::length(cartesian), 0, Math::Epsilon12)) { + return std::optional(); + } + + const double positionX = cartesian.x; + const double positionY = cartesian.y; + const double positionZ = cartesian.z; + + const double oneOverRadiiSquaredX = this->_oneOverRadiiSquared.x; + const double oneOverRadiiSquaredY = this->_oneOverRadiiSquared.y; + const double oneOverRadiiSquaredZ = this->_oneOverRadiiSquared.z; + + const double beta = 1.0 / sqrt( + positionX * positionX * oneOverRadiiSquaredX + + positionY * positionY * oneOverRadiiSquaredY + + positionZ * positionZ * oneOverRadiiSquaredZ); + + return glm::dvec3(positionX * beta, positionY * beta, positionZ * beta); +} + } // namespace CesiumGeospatial diff --git a/CesiumGeospatial/src/SimplePlanarEllipsoidCurve.cpp b/CesiumGeospatial/src/SimplePlanarEllipsoidCurve.cpp new file mode 100644 index 000000000..73cb3337e --- /dev/null +++ b/CesiumGeospatial/src/SimplePlanarEllipsoidCurve.cpp @@ -0,0 +1,106 @@ +#include +#include +#include +#include + +#include + +namespace CesiumGeospatial { + +std::optional +SimplePlanarEllipsoidCurve::fromEarthCenteredEarthFixedCoordinates( + const Ellipsoid& ellipsoid, + const glm::dvec3& sourceEcef, + const glm::dvec3& destinationEcef) { + std::optional scaledSourceEcef = + ellipsoid.scaleToGeocentricSurface(sourceEcef); + std::optional scaledDestinationEcef = + ellipsoid.scaleToGeocentricSurface(destinationEcef); + + if (!scaledSourceEcef.has_value() || !scaledDestinationEcef.has_value()) { + // Unable to scale to geocentric surface coordinates - no curve we can + // generate + return std::optional(); + } + + return SimplePlanarEllipsoidCurve( + ellipsoid, + scaledSourceEcef.value(), + scaledDestinationEcef.value(), + sourceEcef, + destinationEcef); +} + +std::optional +SimplePlanarEllipsoidCurve::fromLongitudeLatitudeHeight( + const Ellipsoid& ellipsoid, + const Cartographic& source, + const Cartographic& destination) { + return SimplePlanarEllipsoidCurve::fromEarthCenteredEarthFixedCoordinates( + ellipsoid, + ellipsoid.cartographicToCartesian(source), + ellipsoid.cartographicToCartesian(destination)); +} + +glm::dvec3 SimplePlanarEllipsoidCurve::getPosition( + double percentage, + double additionalHeight) const { + if (percentage <= 0.0) { + return this->_sourceEcef; + } else if (percentage >= 1.0) { + // We can shortcut our math here and just return the destination. + return this->_destinationEcef; + } + + percentage = glm::clamp(percentage, 0.0, 1.0); + + // Rotate us around the circle between points A and B by the given percentage + // of the total angle we're rotating by. + glm::dvec3 rotatedDirection = + glm::angleAxis(percentage * this->_totalAngle, this->_rotationAxis) * + this->_sourceDirection; + + // It's safe for us to assume here that scaleToGeocentricSurface will return a + // value, since rotatedDirection should never be (0, 0, 0) + glm::dvec3 geocentricPosition = + this->_ellipsoid.scaleToGeocentricSurface(rotatedDirection) + .value_or(glm::dvec3(0, 0, 0)); + + glm::dvec3 geocentricUp = glm::normalize(geocentricPosition); + + double altitudeOffset = + glm::mix(this->_sourceHeight, this->_destinationHeight, percentage) + + additionalHeight; + + return geocentricPosition + geocentricUp * altitudeOffset; +} + +SimplePlanarEllipsoidCurve::SimplePlanarEllipsoidCurve( + const Ellipsoid& ellipsoid, + const glm::dvec3& scaledSourceEcef, + const glm::dvec3& scaledDestinationEcef, + const glm::dvec3& originalSourceEcef, + const glm::dvec3& originalDestinationEcef) + : _ellipsoid(ellipsoid), + _sourceEcef(originalSourceEcef), + _destinationEcef(originalDestinationEcef) { + // Here we find the center of a circle that passes through both the source and + // destination points, and then calculate the angle that we need to move along + // that circle to get from point A to B. + + glm::dquat flyQuat = glm::rotation( + glm::normalize(scaledSourceEcef), + glm::normalize(scaledDestinationEcef)); + + this->_rotationAxis = glm::axis(flyQuat); + this->_totalAngle = glm::angle(flyQuat); + + this->_sourceHeight = glm::length(originalSourceEcef - scaledSourceEcef); + this->_destinationHeight = + glm::length(originalDestinationEcef - scaledDestinationEcef); + + this->_sourceDirection = + glm::normalize(originalSourceEcef - scaledSourceEcef); +} + +} // namespace CesiumGeospatial diff --git a/CesiumGeospatial/test/TestSimplePlanarEllipsoidCurve.cpp b/CesiumGeospatial/test/TestSimplePlanarEllipsoidCurve.cpp new file mode 100644 index 000000000..70ea5f604 --- /dev/null +++ b/CesiumGeospatial/test/TestSimplePlanarEllipsoidCurve.cpp @@ -0,0 +1,213 @@ +#include + +#include +#include + +using namespace CesiumGeospatial; +using namespace CesiumUtility; + +const glm::dvec3 + philadelphiaEcef(1253264.69280105, -4732469.91065521, 4075112.40412297); +const glm::dvec3 + tokyoEcef(-3960158.65587452, 3352568.87555906, 3697235.23506459); +// The antipodal position from the philadelphia coordinates above +const glm::dvec3 philadelphiaAntipodeEcef = + glm::dvec3(-1253369.920224856, 4732412.7444064, -4075146.2160252854); + +// Equivalent to philadelphiaEcef in Longitude, Latitude, Height +const Cartographic philadelphiaLlh( + -1.3119164210487293, + 0.6974930673711344, + 373.64791900173714); +// Equivalent to tokyoEcef +const Cartographic + tokyoLlh(2.4390907007049445, 0.6222806863437318, 283.242432000711); + +TEST_CASE("SimplePlanarEllipsoidCurve::getPosition") { + SECTION("positions at start and end of curve are identical to input " + "coordinates") { + std::optional curve = + SimplePlanarEllipsoidCurve::fromEarthCenteredEarthFixedCoordinates( + Ellipsoid::WGS84, + philadelphiaEcef, + tokyoEcef); + + CHECK(curve.has_value()); + CHECK(Math::equalsEpsilon( + curve.value().getPosition(0.0), + philadelphiaEcef, + Math::Epsilon6)); + CHECK(Math::equalsEpsilon( + curve.value().getPosition(1.0), + tokyoEcef, + Math::Epsilon6)); + } + + SECTION("all points should be coplanar") { + std::optional curve = + SimplePlanarEllipsoidCurve::fromEarthCenteredEarthFixedCoordinates( + Ellipsoid::WGS84, + philadelphiaEcef, + tokyoEcef); + + CHECK(curve.has_value()); + // needs three points to form a plane - get midpoint to make third point + glm::dvec3 midpoint = curve.value().getPosition(0.5); + glm::dvec3 planeNormal = glm::normalize( + glm::cross(philadelphiaEcef - midpoint, tokyoEcef - midpoint)); + + int steps = 100; + + for (int i = 0; i <= steps; i++) { + double time = 1.0 / (double)steps; + double dot = + glm::abs(glm::dot(curve.value().getPosition(time), planeNormal)); + // If the dot product of the point and normal are 0, they're coplanar + CHECK(Math::equalsEpsilon(dot, 0, Math::Epsilon5)); + } + } + + SECTION("should always stay above the earth") { + std::optional curve = + SimplePlanarEllipsoidCurve::fromEarthCenteredEarthFixedCoordinates( + Ellipsoid::WGS84, + philadelphiaEcef, + philadelphiaAntipodeEcef); + + CHECK(curve.has_value()); + + int steps = 100; + for (int i = 0; i <= steps; i++) { + double time = 1.0 / (double)steps; + glm::dvec3 actualResult = curve.value().getPosition(time); + + Cartographic cartographicResult = + Ellipsoid::WGS84.cartesianToCartographic(actualResult) + .value_or(Cartographic(0, 0, 0)); + + CHECK(cartographicResult.height > 0); + } + } + + SECTION("midpoint of reverse path should be identical") { + std::optional forwardCurve = + SimplePlanarEllipsoidCurve::fromEarthCenteredEarthFixedCoordinates( + Ellipsoid::WGS84, + philadelphiaEcef, + tokyoEcef); + + CHECK(forwardCurve.has_value()); + glm::dvec3 forwardResult = forwardCurve.value().getPosition(0.5); + + std::optional reverseCurve = + SimplePlanarEllipsoidCurve::fromEarthCenteredEarthFixedCoordinates( + Ellipsoid::WGS84, + tokyoEcef, + philadelphiaEcef); + + CHECK(reverseCurve.has_value()); + glm::dvec3 reverseResult = reverseCurve.value().getPosition(0.5); + + CHECK(Math::equalsEpsilon(forwardResult, reverseResult, Math::Epsilon6)); + } + + SECTION("curve with same start and end point shouldn't change positions") { + std::optional curve = + SimplePlanarEllipsoidCurve::fromEarthCenteredEarthFixedCoordinates( + Ellipsoid::WGS84, + philadelphiaEcef, + philadelphiaEcef); + + CHECK(curve.has_value()); + + // check a whole bunch of points along the curve to make sure it stays the + // same + const int steps = 25; + for (int i = 0; i <= steps; i++) { + glm::dvec3 result = curve.value().getPosition(i / (double)steps); + CHECK(Math::equalsEpsilon(result, philadelphiaEcef, Math::Epsilon6)); + } + } + + SECTION("should correctly interpolate height") { + double startHeight = 100.0; + double endHeight = 25.0; + + std::optional flightPath = + SimplePlanarEllipsoidCurve::fromLongitudeLatitudeHeight( + Ellipsoid::WGS84, + Cartographic(25.0, 100.0, startHeight), + Cartographic(25.0, 100.0, endHeight)); + + CHECK(flightPath.has_value()); + + std::optional position25Percent = + Ellipsoid::WGS84.cartesianToCartographic( + flightPath.value().getPosition(0.25)); + std::optional position50Percent = + Ellipsoid::WGS84.cartesianToCartographic( + flightPath.value().getPosition(0.5)); + std::optional position75Percent = + Ellipsoid::WGS84.cartesianToCartographic( + flightPath.value().getPosition(0.75)); + + CHECK(position25Percent.has_value()); + CHECK(Math::equalsEpsilon( + position25Percent.value().height, + (endHeight - startHeight) * 0.25 + startHeight, + Math::Epsilon6)); + CHECK(position50Percent.has_value()); + CHECK(Math::equalsEpsilon( + position50Percent.value().height, + (endHeight - startHeight) * 0.5 + startHeight, + Math::Epsilon6)); + CHECK(position75Percent.has_value()); + CHECK(Math::equalsEpsilon( + position75Percent.value().height, + (endHeight - startHeight) * 0.75 + startHeight, + Math::Epsilon6)); + } +} + +TEST_CASE( + "SimplePlanarEllipsoidCurve::fromEarthCenteredEarthFixedCoordinates") { + SECTION("should fail on coordinates (0, 0, 0)") { + std::optional curve = + SimplePlanarEllipsoidCurve::fromEarthCenteredEarthFixedCoordinates( + Ellipsoid::WGS84, + philadelphiaEcef, + glm::dvec3(0, 0, 0)); + + CHECK(!curve.has_value()); + } +} + +TEST_CASE("SimplePlanarEllipsoidCurve::fromLongitudeLatitudeHeight") { + SECTION("should match results of curve created from equivalent ECEF " + "coordinates") { + std::optional llhCurve = + SimplePlanarEllipsoidCurve::fromLongitudeLatitudeHeight( + Ellipsoid::WGS84, + philadelphiaLlh, + tokyoLlh); + + std::optional ecefCurve = + SimplePlanarEllipsoidCurve::fromEarthCenteredEarthFixedCoordinates( + Ellipsoid::WGS84, + philadelphiaEcef, + tokyoEcef); + + CHECK(llhCurve.has_value()); + CHECK(ecefCurve.has_value()); + + int steps = 100; + + for (int i = 0; i <= steps; i++) { + double n = 1.0 / (double)steps; + CHECK(Math::equalsEpsilon( + ecefCurve.value().getPosition(n), + llhCurve.value().getPosition(n), + Math::Epsilon6)); + } + } +}