diff --git a/CHANGES.md b/CHANGES.md index c54d72428..f0b62d2d3 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -6,6 +6,8 @@ - Added `convertAccessorTypeToPropertyType` and `convertPropertyTypeToAccessorType` to `CesiumGltf::PropertyType`. - Added support for building in `vcpkg` manifest mode. +- Added support for orthographic and skewed perspective views. +- Added an overload of `Math::equalsEpsilon` for glm matrices. ##### Fixes :wrench: diff --git a/Cesium3DTilesSelection/include/Cesium3DTilesSelection/ViewState.h b/Cesium3DTilesSelection/include/Cesium3DTilesSelection/ViewState.h index 1ff442a03..7c5ae0d84 100644 --- a/Cesium3DTilesSelection/include/Cesium3DTilesSelection/ViewState.h +++ b/Cesium3DTilesSelection/include/Cesium3DTilesSelection/ViewState.h @@ -8,9 +8,11 @@ #include #include +#include #include #include +#include #include namespace Cesium3DTilesSelection { @@ -22,10 +24,23 @@ namespace Cesium3DTilesSelection { */ class CESIUM3DTILESSELECTION_API ViewState final { - // TODO: Add support for orthographic and off-center perspective frustums public: /** - * @brief Creates a new instance of a view state. + * @brief Creates a new instance of a view state with a symmetric perspective + * projection. + * @deprecated Use {@link ViewState::ViewState} instead. + */ + [[deprecated("Use ViewState::ViewState instead.")]] static ViewState create( + const glm::dvec3& position, + const glm::dvec3& direction, + const glm::dvec3& up, + const glm::dvec2& viewportSize, + double horizontalFieldOfView, + double verticalFieldOfView, + const CesiumGeospatial::Ellipsoid& ellipsoid CESIUM_DEFAULT_ELLIPSOID); + + /** + * @brief Creates a new instance. * * @param position The position of the eye point of the camera. * @param direction The view direction vector of the camera. @@ -40,7 +55,7 @@ class CESIUM3DTILESSELECTION_API ViewState final { * from the cartesian position. * Default value: {@link CesiumGeospatial::Ellipsoid::WGS84}. */ - static ViewState create( + ViewState( const glm::dvec3& position, const glm::dvec3& direction, const glm::dvec3& up, @@ -49,6 +64,53 @@ class CESIUM3DTILESSELECTION_API ViewState final { double verticalFieldOfView, const CesiumGeospatial::Ellipsoid& ellipsoid CESIUM_DEFAULT_ELLIPSOID); + /** + * @brief Creates a new instance of a view state with a general projection, + * including skewed perspective and orthographic projections. + * + * @param viewMatrix The view's view matrix i.e., the inverse of its pose + * @param projectionMatrix see e.g. + * {@link CesiumGeometry::Transforms::createPerspectiveMatrix} + * @param viewportSize The size of the viewport, in pixels. + * @param ellipsoid The ellipsoid that will be used to compute the + * {@link ViewState#getPositionCartographic cartographic position} + * from the cartesian position. + * Default value: {@link CesiumGeospatial::Ellipsoid::WGS84}. + */ + ViewState( + const glm::dmat4& viewMatrix, + const glm::dmat4& projectionMatrix, + const glm::dvec2& viewportSize, + const CesiumGeospatial::Ellipsoid& ellipsoid CESIUM_DEFAULT_ELLIPSOID); + + /** + * @brief Creates a new instance of a view state with an orthographic + * projection. + * + * @param position The position of the eye point of the camera. + * @param direction The view direction vector of the camera. + * @param up The up vector of the camera. + * @param viewportSize The size of the viewport, in pixels. + * @param left left distance of near plane edge from center + * @param right right distance of near plane edge + * @param bottom bottom distance of near plane edge + * @param top top distance of near plane edge + * @param ellipsoid The ellipsoid that will be used to compute the + * {@link ViewState#getPositionCartographic cartographic position} + * from the cartesian position. + * Default value: {@link CesiumGeospatial::Ellipsoid::WGS84}. + */ + ViewState( + const glm::dvec3& position, + const glm::dvec3& direction, + const glm::dvec3& up, + const glm::dvec2& viewportSize, + double left, + double right, + double bottom, + double top, + const CesiumGeospatial::Ellipsoid& ellipsoid CESIUM_DEFAULT_ELLIPSOID); + /** * @brief Gets the position of the camera in Earth-centered, Earth-fixed * coordinates. @@ -65,7 +127,12 @@ class CESIUM3DTILESSELECTION_API ViewState final { * @brief Gets the up direction of the camera in Earth-centered, Earth-fixed * coordinates. */ - const glm::dvec3& getUp() const noexcept { return this->_up; } + glm::dvec3 getUp() const noexcept { + return { + this->_viewMatrix[0][1], + this->_viewMatrix[1][1], + this->_viewMatrix[2][1]}; + } /** * @brief Gets the position of the camera as a longitude / latitude / height. @@ -86,17 +153,27 @@ class CESIUM3DTILESSELECTION_API ViewState final { } /** - * @brief Gets the horizontal field-of-view angle in radians. + * @brief Gets the horizontal field-of-view angle in radians. This is only + * valid for a symmetric perspective projection. */ - double getHorizontalFieldOfView() const noexcept { - return this->_horizontalFieldOfView; - } + double getHorizontalFieldOfView() const noexcept; /** - * @brief Gets the vertical field-of-view angle in radians. + * @brief Gets the vertical field-of-view angle in radians. This is only + * valid for a symmetric perspective projection. */ - double getVerticalFieldOfView() const noexcept { - return this->_verticalFieldOfView; + double getVerticalFieldOfView() const noexcept; + + /** + * @brief Gets the view matrix for the ViewState. + */ + const glm::dmat4& getViewMatrix() const noexcept { return this->_viewMatrix; } + + /** + * @brief Gets the projection matrix for the ViewState. + */ + const glm::dmat4& getProjectionMatrix() const noexcept { + return this->_projectionMatrix; } /** @@ -142,40 +219,16 @@ class CESIUM3DTILESSELECTION_API ViewState final { const noexcept; private: - /** - * @brief Creates a new instance. - * - * @param position The position of the eye point of the camera. - * @param direction The view direction vector of the camera. - * @param up The up vector of the camera. - * @param viewportSize The size of the viewport, in pixels. - * @param horizontalFieldOfView The horizontal field-of-view (opening) - * angle of the camera, in radians. - * @param verticalFieldOfView The vertical field-of-view (opening) - * angle of the camera, in radians. - */ - ViewState( - const glm::dvec3& position, - const glm::dvec3& direction, - const glm::dvec3& up, - const glm::dvec2& viewportSize, - double horizontalFieldOfView, - double verticalFieldOfView, - const std::optional& positionCartographic, - const CesiumGeospatial::Ellipsoid& ellipsoid); - const glm::dvec3 _position; const glm::dvec3 _direction; - const glm::dvec3 _up; const glm::dvec2 _viewportSize; - const double _horizontalFieldOfView; - const double _verticalFieldOfView; const CesiumGeospatial::Ellipsoid _ellipsoid; - const double _sseDenominator; const std::optional _positionCartographic; const CesiumGeometry::CullingVolume _cullingVolume; + const glm::dmat4 _viewMatrix; + const glm::dmat4 _projectionMatrix; }; } // namespace Cesium3DTilesSelection diff --git a/Cesium3DTilesSelection/src/ViewState.cpp b/Cesium3DTilesSelection/src/ViewState.cpp index 11273bb10..ea3ed3bb2 100644 --- a/Cesium3DTilesSelection/src/ViewState.cpp +++ b/Cesium3DTilesSelection/src/ViewState.cpp @@ -5,17 +5,21 @@ #include #include #include +#include #include #include -#include #include #include #include +#include #include #include -#include +#include +#include +#include +#include #include #include @@ -39,10 +43,27 @@ namespace Cesium3DTilesSelection { viewportSize, horizontalFieldOfView, verticalFieldOfView, - ellipsoid.cartesianToCartographic(position), ellipsoid); } +namespace { +glm::dvec3 positionFromView(const glm::dmat4& viewMatrix) { + // Back out the world position by multiplying the view matrix translation by + // the rotation transpose (inverse) and negating. + glm::dvec3 position(0.0); + position.x = -glm::dot(glm::dvec3(viewMatrix[0]), glm::dvec3(viewMatrix[3])); + position.y = -glm::dot(glm::dvec3(viewMatrix[1]), glm::dvec3(viewMatrix[3])); + position.z = -glm::dot(glm::dvec3(viewMatrix[2]), glm::dvec3(viewMatrix[3])); + return position; +} + +glm::dvec3 directionFromView(const glm::dmat4& viewMatrix) { + return glm::dvec3(viewMatrix[0][2], viewMatrix[1][2], viewMatrix[2][2]) * + -1.0; +} +} // namespace + +// The near and far plane values don't matter for frustum culling. ViewState::ViewState( const glm::dvec3& position, const glm::dvec3& direction, @@ -50,23 +71,53 @@ ViewState::ViewState( const glm::dvec2& viewportSize, double horizontalFieldOfView, double verticalFieldOfView, - const std::optional& positionCartographic, const CesiumGeospatial::Ellipsoid& ellipsoid) - : _position(position), - _direction(direction), - _up(up), + : ViewState( + Transforms::createViewMatrix(position, direction, up), + Transforms::createPerspectiveMatrix( + horizontalFieldOfView, + verticalFieldOfView, + 100.0, + std::numeric_limits::infinity()), + viewportSize, + ellipsoid) {} + +ViewState::ViewState( + const glm::dmat4& viewMatrix, + const glm::dmat4& projectionMatrix, + const glm::dvec2& viewportSize, + const CesiumGeospatial::Ellipsoid& ellipsoid) + : _position(positionFromView(viewMatrix)), + _direction(directionFromView(viewMatrix)), _viewportSize(viewportSize), - _horizontalFieldOfView(horizontalFieldOfView), - _verticalFieldOfView(verticalFieldOfView), _ellipsoid(ellipsoid), - _sseDenominator(2.0 * glm::tan(0.5 * verticalFieldOfView)), - _positionCartographic(positionCartographic), - _cullingVolume(createCullingVolume( - position, - direction, - up, - horizontalFieldOfView, - verticalFieldOfView)) {} + _positionCartographic( + ellipsoid.cartesianToCartographic(positionFromView(viewMatrix))), + _cullingVolume(createCullingVolume(projectionMatrix * viewMatrix)), + _viewMatrix(viewMatrix), + _projectionMatrix(projectionMatrix) {} + +ViewState::ViewState( + const glm::dvec3& position, + const glm::dvec3& direction, + const glm::dvec3& up, + const glm::dvec2& viewportSize, + double left, + double right, + double bottom, + double top, + const CesiumGeospatial::Ellipsoid& ellipsoid) + : ViewState( + Transforms::createViewMatrix(position, direction, up), + Transforms::createOrthographicMatrix( + left, + right, + bottom, + top, + 100.0, + std::numeric_limits::infinity()), + viewportSize, + ellipsoid) {} namespace { template @@ -205,7 +256,29 @@ double ViewState::computeScreenSpaceError( double distance) const noexcept { // Avoid divide by zero when viewer is inside the tile distance = glm::max(distance, 1e-7); - const double sseDenominator = this->_sseDenominator; - return (geometricError * this->_viewportSize.y) / (distance * sseDenominator); + // This is a simplified version of the projection transform and homogeneous + // division. We transform the coordinate (0.0, geometric error, -distance, + // 1) and use the resulting NDC to find the screen space error. That's not + // quite right: the distance argument is actually the slant distance to the + // tile, whereas view-space Z is perpendicular to the near plane. + const glm::dmat4& projMat = this->_projectionMatrix; + glm::dvec4 centerNdc = projMat * glm::dvec4(0.0, 0.0, -distance, 1.0); + centerNdc /= centerNdc.w; + glm::dvec4 errorOffsetNdc = + projMat * glm::dvec4(0.0, geometricError, -distance, 1.0); + errorOffsetNdc /= errorOffsetNdc.w; + + double ndcError = (errorOffsetNdc - centerNdc).y; + // ndc bounds are [-1.0, 1.0]. Our projection matrix has the top of the + // screen at -1.0, so ndcError is negative. + return -ndcError * this->_viewportSize.y / 2.0; +} + +double ViewState::getHorizontalFieldOfView() const noexcept { + return std::atan(1.0 / this->_projectionMatrix[0][0]) * 2.0; +} + +double ViewState::getVerticalFieldOfView() const noexcept { + return std::atan(-1.0 / this->_projectionMatrix[1][1]) * 2.0; } } // namespace Cesium3DTilesSelection diff --git a/Cesium3DTilesSelection/test/TestTilesetSelectionAlgorithm.cpp b/Cesium3DTilesSelection/test/TestTilesetSelectionAlgorithm.cpp index 18a0bbfb9..ee5213242 100644 --- a/Cesium3DTilesSelection/test/TestTilesetSelectionAlgorithm.cpp +++ b/Cesium3DTilesSelection/test/TestTilesetSelectionAlgorithm.cpp @@ -85,7 +85,7 @@ static void initializeTileset(Tileset& tileset) { double horizontalFieldOfView = Math::degreesToRadians(60.0); double verticalFieldOfView = std::atan(std::tan(horizontalFieldOfView * 0.5) / aspectRatio) * 2.0; - ViewState viewState = ViewState::create( + ViewState viewState = ViewState( viewPosition, glm::normalize(viewFocus - viewPosition), viewUp, @@ -117,7 +117,7 @@ static ViewState zoomToTile(const Tile& tile) { double horizontalFieldOfView = Math::degreesToRadians(60.0); double verticalFieldOfView = std::atan(std::tan(horizontalFieldOfView * 0.5) / aspectRatio) * 2.0; - return ViewState::create( + return ViewState( viewPosition, glm::normalize(viewFocus - viewPosition), viewUp, @@ -206,7 +206,7 @@ TEST_CASE("Test replace refinement for render") { // Zoom out from the tileset a little bit to make sure the root meet sse glm::dvec3 zoomOutPosition = viewState.getPosition() - viewState.getDirection() * 2500.0; - ViewState zoomOutViewState = ViewState::create( + ViewState zoomOutViewState = ViewState( zoomOutPosition, viewState.getDirection(), viewState.getUp(), @@ -330,7 +330,7 @@ TEST_CASE("Test replace refinement for render") { ViewState viewState = zoomToTileset(tileset); glm::dvec3 zoomInPosition = viewState.getPosition() + viewState.getDirection() * 200.0; - ViewState zoomInViewState = ViewState::create( + ViewState zoomInViewState = ViewState( zoomInPosition, viewState.getDirection(), viewState.getUp(), @@ -422,7 +422,7 @@ TEST_CASE("Test replace refinement for render") { { glm::dvec3 zoomOutPosition = viewState.getPosition() - viewState.getDirection() * 100.0; - ViewState zoomOutViewState = ViewState::create( + ViewState zoomOutViewState = ViewState( zoomOutPosition, viewState.getDirection(), viewState.getUp(), @@ -850,7 +850,7 @@ TEST_CASE("Test multiple frustums") { // Zoom out from the tileset a little bit to make sure the root meets sse glm::dvec3 zoomOutPosition = viewState.getPosition() - viewState.getDirection() * 2500.0; - ViewState zoomOutViewState = ViewState::create( + ViewState zoomOutViewState = ViewState( zoomOutPosition, viewState.getDirection(), viewState.getUp(), @@ -922,7 +922,7 @@ TEST_CASE("Test multiple frustums") { // (child of the first child of the root). glm::dvec3 zoomInPosition = zoomToTileViewState.getPosition() + zoomToTileViewState.getDirection() * 250.0; - ViewState zoomInViewState1 = ViewState::create( + ViewState zoomInViewState1 = ViewState( zoomInPosition, zoomToTileViewState.getDirection(), zoomToTileViewState.getUp(), @@ -934,7 +934,7 @@ TEST_CASE("Test multiple frustums") { zoomInPosition = zoomToTileViewState.getPosition() + glm::dvec3(15.0, 0, 0) + zoomToTileViewState.getDirection() * 243.0; - ViewState zoomInViewState2 = ViewState::create( + ViewState zoomInViewState2 = ViewState( zoomInPosition, zoomToTileViewState.getDirection(), zoomToTileViewState.getUp(), @@ -1705,7 +1705,7 @@ TEST_CASE("Additive-refined tiles are added to the tilesFadingOut array") { REQUIRE(position); position->height += 100000; - ViewState zoomedOut = ViewState::create( + ViewState zoomedOut = ViewState( Ellipsoid::WGS84.cartographicToCartesian(*position), viewState.getDirection(), viewState.getUp(), diff --git a/CesiumGeometry/include/CesiumGeometry/CullingVolume.h b/CesiumGeometry/include/CesiumGeometry/CullingVolume.h index c5ce5157a..7434931ce 100644 --- a/CesiumGeometry/include/CesiumGeometry/CullingVolume.h +++ b/CesiumGeometry/include/CesiumGeometry/CullingVolume.h @@ -2,6 +2,8 @@ #include +#include + namespace CesiumGeometry { /** @@ -58,4 +60,61 @@ CullingVolume createCullingVolume( const glm::dvec3& up, double fovx, double fovy) noexcept; + +/** + * @brief Create a {@link CullingVolume} from a projection matrix. + * + * The matrix can be a composite view - projection matrix; the volume will then + * cull world coordinates. It can also be a model - view - projection matrix, + * which will cull model coordinates. + */ +CullingVolume createCullingVolume(const glm::dmat4& clipMatrix); + +/** + * @brief Creates a {@link CullingVolume} for a perspective frustum. + * + * @param position The eye position + * @param direction The viewing direction + * @param up The up-vector of the frustum + * @param l left edge + * @param r right edge + * @param t top edge + * @param b bottom edge + * @param n near plane distance + * @return The {@link CullingVolume} + */ +CullingVolume createCullingVolume( + const glm::dvec3& position, + const glm::dvec3& direction, + const glm::dvec3& up, + double l, + double r, + double b, + double t, + double n) noexcept; + +/** + * @brief Creates a {@link CullingVolume} for an orthographic frustum. + * + * @param position The eye position + * @param direction The viewing direction + * @param up The up-vector of the frustum + * @param l left edge + * @param r right edge + * @param t top edge + * @param b bottom edge + * @param n near plane distance + * @return The {@link CullingVolume} + */ + +CullingVolume createOrthographicCullingVolume( + const glm::dvec3& position, + const glm::dvec3& direction, + const glm::dvec3& up, + double l, + double r, + double b, + double t, + double n) noexcept; + } // namespace CesiumGeometry diff --git a/CesiumGeometry/include/CesiumGeometry/Transforms.h b/CesiumGeometry/include/CesiumGeometry/Transforms.h index 9e1c6fee1..41b8e7cd8 100644 --- a/CesiumGeometry/include/CesiumGeometry/Transforms.h +++ b/CesiumGeometry/include/CesiumGeometry/Transforms.h @@ -95,6 +95,91 @@ struct CESIUMGEOMETRY_API Transforms final { */ static const glm::dmat4& getUpAxisTransform(CesiumGeometry::Axis from, CesiumGeometry::Axis to); + + /** + * @brief Create a view matrix. + * + * This is similar to glm::lookAt(), but uses the pose of the viewer to create + * the view matrix. The view matrix is the inverse of the pose matrix. + * + * @param position position of the eye + * @param direction view vector i.e., -z axis of the viewer's pose. + * @param up up vector of viewer i.e., y axis of the viewer's pose. + * + * @returns The view matrix. + */ + static glm::dmat4 createViewMatrix( + const glm::dvec3& position, + const glm::dvec3& direction, + const glm::dvec3& up); + + /** + * @brief Compute a Vulkan-style perspective projection matrix with reversed + * Z. + * + * "Vulkan-style", as return by this function and others, uses the following + * conventions: + * * X maps from -1 to 1 left to right + * * Y maps from 1 to -1 bottom to top + * * Z maps from 1 to 0 near to far (known as "reverse Z") + * + * @param fovx horizontal field of view in radians + * @param fovy vertical field of view in radians + * @param zNear distance to near plane + * @param zFar distance to far plane + */ + static glm::dmat4 + createPerspectiveMatrix(double fovx, double fovy, double zNear, double zFar); + + /** + * @brief Compute a Vulkan-style perspective projection matrix with reversed + * Z. + * + * "Vulkan-style", as return by this function and others, uses the following + * conventions: + * * X maps from -1 to 1 left to right + * * Y maps from 1 to -1 bottom to top + * * Z maps from 1 to 0 near to far (known as "reverse Z") + * + * @param left left distance of near plane edge from center + * @param right right distance of near plane edge + * @param bottom bottom distance of near plane edge + * @param top top distance of near plane edge + * @param zNear distance of near plane + * @param zFar distance of far plane. This can be infinite + */ + static glm::dmat4 createPerspectiveMatrix( + double left, + double right, + double bottom, + double top, + double zNear, + double zFar); + + /** + * @brief Compute a Vulkan-style orthographic projection matrix with reversed + * Z. + * + * "Vulkan-style", as return by this function and others, uses the following + * conventions: + * * X maps from -1 to 1 left to right + * * Y maps from 1 to -1 bottom to top + * * Z maps from 1 to 0 near to far (known as "reverse Z") + * + * @param left left distance of near plane edge from center + * @param right right distance of near plane edge + * @param bottom bottom distance of near plane edge + * @param top top distance of near plane edge + * @param zNear distance of near plane + * @param zFar distance of far plane. This can be infinite + */ + static glm::dmat4 createOrthographicMatrix( + double left, + double right, + double bottom, + double top, + double zNear, + double zFar); }; } // namespace CesiumGeometry diff --git a/CesiumGeometry/src/CullingVolume.cpp b/CesiumGeometry/src/CullingVolume.cpp index ee11a4f7c..960881437 100644 --- a/CesiumGeometry/src/CullingVolume.cpp +++ b/CesiumGeometry/src/CullingVolume.cpp @@ -1,6 +1,9 @@ #include #include +#include +#include +#include #include #include #include @@ -75,4 +78,89 @@ CullingVolume createCullingVolume( return {leftPlane, rightPlane, topPlane, bottomPlane}; } + +CullingVolume createCullingVolume(const glm::dmat4& clipMatrix) { + // Gribb / Hartmann method. + // See + // https://www.gamedevs.org/uploads/fast-extraction-viewing-frustum-planes-from-world-view-projection-matrix.pdf + // Here's a sketch of the method for the left clipping plane. Given: + // v incoming eye-space vertex + // x', w' coordinate values after matrix multiplication + // + // The clipping test is -w' < x' + // or -(row3 dot v) < (row0 dot v) + // 0 < (row3 dot v) + (row0 dot v) + // 0 < (row3 + row0) dot v + // This is a test against the plane a*x + b*y + c*z + d*w = 0, but the + // incoming w = 1, so we have a standard plane a*x + b*y + c*z + d = 0 where + // the factors are the sums of the elements of the matrix rows. + // + auto normalizePlane = [](double a, double b, double c, double d) { + double len = sqrt(a * a + b * b + c * c); + return Plane(glm::dvec3(a / len, b / len, c / len), d / len); + }; + const Plane leftPlane = normalizePlane( + clipMatrix[0][3] + clipMatrix[0][0], + clipMatrix[1][3] + clipMatrix[1][0], + clipMatrix[2][3] + clipMatrix[2][0], + clipMatrix[3][3] + clipMatrix[3][0]); + const Plane rightPlane = normalizePlane( + clipMatrix[0][3] - clipMatrix[0][0], + clipMatrix[1][3] - clipMatrix[1][0], + clipMatrix[2][3] - clipMatrix[2][0], + clipMatrix[3][3] - clipMatrix[3][0]); + const Plane bottomPlane = normalizePlane( + clipMatrix[0][3] - clipMatrix[0][1], + clipMatrix[1][3] - clipMatrix[1][1], + clipMatrix[2][3] - clipMatrix[2][1], + clipMatrix[3][3] - clipMatrix[3][1]); + const Plane topPlane = normalizePlane( + clipMatrix[0][3] + clipMatrix[0][1], + clipMatrix[1][3] + clipMatrix[1][1], + clipMatrix[2][3] + clipMatrix[2][1], + clipMatrix[3][3] + clipMatrix[3][1]); + return {leftPlane, rightPlane, topPlane, bottomPlane}; +} + +CullingVolume createCullingVolume( + const glm::dvec3& position, + const glm::dvec3& direction, + const glm::dvec3& up, + double l, + double r, + double b, + double t, + double n) noexcept { + glm::dmat4 projMatrix = Transforms::createPerspectiveMatrix( + l, + r, + b, + t, + n, + std::numeric_limits::infinity()); + glm::dmat4 viewMatrix = glm::lookAt(position, position + direction, up); + glm::dmat4 clipMatrix = projMatrix * viewMatrix; + return createCullingVolume(clipMatrix); +} + +CullingVolume createOrthographicCullingVolume( + const glm::dvec3& position, + const glm::dvec3& direction, + const glm::dvec3& up, + double l, + double r, + double b, + double t, + double n) noexcept { + glm::dmat4 projMatrix = Transforms::createOrthographicMatrix( + l, + r, + b, + t, + n, + std::numeric_limits::infinity()); + glm::dmat4 viewMatrix = glm::lookAt(position, position + direction, up); + glm::dmat4 clipMatrix = projMatrix * viewMatrix; + return createCullingVolume(clipMatrix); +} } // namespace CesiumGeometry diff --git a/CesiumGeometry/src/Transforms.cpp b/CesiumGeometry/src/Transforms.cpp index 3a92dd9e0..6f9fe16a4 100644 --- a/CesiumGeometry/src/Transforms.cpp +++ b/CesiumGeometry/src/Transforms.cpp @@ -5,6 +5,9 @@ #include #include +#include +#include + namespace CesiumGeometry { namespace { @@ -157,4 +160,104 @@ const glm::dmat4& Transforms::getUpAxisTransform(Axis from, Axis to) { return Identity; } +glm::dmat4 Transforms::createViewMatrix( + const glm::dvec3& position, + const glm::dvec3& direction, + const glm::dvec3& up) { + const glm::dvec3 forward = -1.0 * direction; + const glm::dvec3 side = glm::normalize(glm::cross(up, forward)); + const glm::dvec3 poseUp = glm::normalize(glm::cross(forward, side)); + // Invert the pose to create the view matrix. Transpose the rotation part and + // invert the position "by hand." + glm::dmat4 result(1.0); + result[0][0] = side.x; + result[1][0] = side.y; + result[2][0] = side.z; + result[0][1] = poseUp.x; + result[1][1] = poseUp.y; + result[2][1] = poseUp.z; + result[0][2] = forward.x; + result[1][2] = forward.y; + result[2][2] = forward.z; + result[3][0] = -glm::dot(side, position); + result[3][1] = -glm::dot(poseUp, position); + result[3][2] = -glm::dot(forward, position); + return result; +} + +glm::dmat4 Transforms::createPerspectiveMatrix( + double fovx, + double fovy, + double zNear, + double zFar) { + double m22; + double m32; + if (zFar == std::numeric_limits::infinity()) { + m22 = 0.0; + m32 = zNear; + } else { + m22 = zNear / (zFar - zNear); + m32 = zNear * zFar / (zFar - zNear); + } + return glm::dmat4( + glm::dvec4(1.0 / std::tan(fovx * 0.5), 0.0, 0.0, 0.0), + glm::dvec4(0.0, -1.0 / std::tan(fovy * 0.5), 0.0, 0.0), + glm::dvec4(0.0, 0.0, m22, -1.0), + glm::dvec4(0.0, 0.0, m32, 0.0)); +} + +glm::dmat4 Transforms::createPerspectiveMatrix( + double left, + double right, + double bottom, + double top, + double zNear, + double zFar) { + double m22; + double m32; + if (zFar == std::numeric_limits::infinity()) { + m22 = 0.0; + m32 = zNear; + } else { + m22 = zNear / (zFar - zNear); + m32 = zNear * zFar / (zFar - zNear); + } + return glm::dmat4( + glm::dvec4(2.0 * zNear / (right - left), 0.0, 0.0, 0.0), + glm::dvec4(0.0, 2.0 * zNear / (bottom - top), 0.0, 0.0), + glm::dvec4( + (right + left) / (right - left), + (bottom + top) / (bottom - top), + m22, + -1.0), + glm::dvec4(0.0, 0.0, m32, 0.0)); +} + +glm::dmat4 Transforms::createOrthographicMatrix( + double left, + double right, + double bottom, + double top, + double zNear, + double zFar) { + double m22; + double m32; + if (zFar == std::numeric_limits::infinity()) { + m22 = 0.0; + m32 = 1.0; + } else { + m22 = 1.0 / (zFar - zNear); + m32 = zFar / (zFar - zNear); + } + return glm::dmat4( + glm::dvec4(2.0 / (right - left), 0.0, 0.0, 0.0), + glm::dvec4(0.0, 2.0 / (bottom - top), 0.0, 0.0), + glm::dvec4(0.0, 0.0, m22, 0.0), + glm::dvec4( + -(right + left) / (right - left), + -(bottom + top) / (bottom - top), + m32, + 1.0)); +} + } // namespace CesiumGeometry diff --git a/CesiumGeometry/test/TestCullingVolume.cpp b/CesiumGeometry/test/TestCullingVolume.cpp index a6c470f7e..5a6fa8652 100644 --- a/CesiumGeometry/test/TestCullingVolume.cpp +++ b/CesiumGeometry/test/TestCullingVolume.cpp @@ -1,4 +1,5 @@ #include +#include #include #include @@ -13,6 +14,32 @@ using namespace doctest; using namespace CesiumGeometry; using namespace CesiumUtility; +namespace { +bool equalsEpsilon( + const Plane& left, + const Plane& right, + double relativeEpsilon) { + return Math::equalsEpsilon( + left.getNormal(), + right.getNormal(), + relativeEpsilon) && + Math::equalsEpsilon( + left.getDistance(), + right.getDistance(), + relativeEpsilon); +} +} // namespace + +bool equalsEpsilon( + const CullingVolume& left, + const CullingVolume& right, + double relativeEpsilon) { + return equalsEpsilon(left.leftPlane, right.leftPlane, relativeEpsilon) && + equalsEpsilon(left.rightPlane, right.rightPlane, relativeEpsilon) && + equalsEpsilon(left.topPlane, right.topPlane, relativeEpsilon) && + equalsEpsilon(left.bottomPlane, right.bottomPlane, relativeEpsilon); +} + TEST_CASE("CullingVolume::createCullingVolume") { SUBCASE("Shouldn't crash when too far from the globe") { CHECK_NOTHROW(createCullingVolume( @@ -31,4 +58,26 @@ TEST_CASE("CullingVolume::createCullingVolume") { Math::PiOverTwo, Math::PiOverTwo)); } -} \ No newline at end of file +} + +TEST_CASE("CullingVolume construction") { + SUBCASE("Field of view and clip matrix") { + glm::dvec3 position(1e5, 1e5, 1e5); + glm::dvec3 direction(0, 0, 1); + glm::dvec3 up(0, 1, 0); + CullingVolume traditional = createCullingVolume( + position, + direction, + up, + Math::PiOverTwo, + Math::PiOverTwo); + CullingVolume rad = createCullingVolume( + Transforms::createPerspectiveMatrix( + Math::PiOverTwo, + Math::PiOverTwo, + 10, + 200000) * + Transforms::createViewMatrix(position, direction, up)); + CHECK(equalsEpsilon(traditional, rad, 1e-10)); + } +} diff --git a/CesiumGeometry/test/TestTransforms.cpp b/CesiumGeometry/test/TestTransforms.cpp index b11cc05fb..bae29027e 100644 --- a/CesiumGeometry/test/TestTransforms.cpp +++ b/CesiumGeometry/test/TestTransforms.cpp @@ -1,10 +1,15 @@ #include #include +#include #include #include #include +#include +#include +#include + using namespace CesiumGeometry; TEST_CASE("Transforms convert the axes correctly") { @@ -98,3 +103,113 @@ TEST_CASE("Gets up axis transform") { CHECK(Transforms::getUpAxisTransform(Axis::Z, Axis::Z) == Identity); } } + +namespace { +bool pointInClipVolume(const glm::dvec4& point) { + const double w = point.w; + return -w <= point.x && point.x <= w && -w <= point.y && point.y <= w && + 0.0 <= point.z && point.z <= w; +} +} // namespace +TEST_CASE("Test perspective projection matrices") { + const double horizontalFieldOfView = + CesiumUtility::Math::degreesToRadians(60.0); + const double verticalFieldOfView = + CesiumUtility::Math::degreesToRadians(45.0); + const double zNear = 1.0; + const double zFar = 20000.0; + const glm::dmat4 projMat = Transforms::createPerspectiveMatrix( + horizontalFieldOfView, + verticalFieldOfView, + zNear, + zFar); + std::vector testPoints; + testPoints.reserve(1000); + for (int i = 0; i < 11; ++i) { + double hRad = + -horizontalFieldOfView / 2.0 + i * horizontalFieldOfView / 10.0; + hRad = CesiumUtility::Math::clamp( + hRad, + -horizontalFieldOfView + 0.1, + horizontalFieldOfView - 0.1); + double sinH = std::sin(hRad); + for (int j = 0; j < 10; ++j) { + double vRad = -verticalFieldOfView / 2.0 + j * verticalFieldOfView / 10.0; + vRad = CesiumUtility::Math::clamp( + vRad, + -verticalFieldOfView + 0.1, + verticalFieldOfView - 0.1); + double sinV = std::sin(vRad); + for (int k = 0; k < 10; ++k) { + double z = zNear + k * (zFar - zNear) / 10.0; + z = CesiumUtility::Math::clamp(z, zNear + .1, zFar - .1); + testPoints.emplace_back(sinH * z, sinV * z, -z, 1.0); + } + } + } + SUBCASE("Check that points lie in clipping volume") { + CHECK(std::all_of( + testPoints.begin(), + testPoints.end(), + [&projMat](const glm::dvec4& p) { + glm::dvec4 projected = projMat * p; + return pointInClipVolume(projected); + })); + } + double hDim = std::tan(horizontalFieldOfView / 2.0) * zNear; + double vDim = std::tan(verticalFieldOfView / 2.0) * zNear; + glm::dmat4 corners = Transforms::createPerspectiveMatrix( + -hDim, + hDim, + -vDim, + vDim, + zNear, + zFar); + SUBCASE("Check that different perspective constructions are equivalent") { + CHECK(CesiumUtility::Math::equalsEpsilon(projMat, corners, 1e-14, 1e-14)); + } + + SUBCASE("Check skewed projection") { + glm::dmat4 skewed = + Transforms::createPerspectiveMatrix(0.0, hDim, 0.0, vDim, zNear, zFar); + for (auto& point : testPoints) { + glm::dvec4 skewProjected = skewed * point; + if (pointInClipVolume(skewProjected)) { + glm::dvec4 symProjected = corners * point; + skewProjected /= skewProjected.w; + symProjected /= symProjected.w; + CHECK(CesiumUtility::Math::equalsEpsilon( + skewProjected.x / 2.0 + .5, + symProjected.x, + 1e-14)); + CHECK(CesiumUtility::Math::equalsEpsilon( + skewProjected.y / 2.0 - .5, + symProjected.y, + 1e-14)); + CHECK(CesiumUtility::Math::equalsEpsilon( + skewProjected.z, + symProjected.z, + 1e-14)); + } + } + } + + SUBCASE("Check that points are contained in orthgraphic projection too") { + double orthoHDim = hDim / zNear * zFar; + double orthoVDim = vDim / zNear * zFar; + glm::dmat4 ortho = Transforms::createOrthographicMatrix( + -orthoHDim, + orthoHDim, + -orthoVDim, + orthoVDim, + zNear, + zFar); + CHECK(std::all_of( + testPoints.begin(), + testPoints.end(), + [&ortho](const glm::dvec4& p) { + glm::dvec4 projected = ortho * p; + return pointInClipVolume(projected); + })); + } +} diff --git a/CesiumUtility/include/CesiumUtility/Math.h b/CesiumUtility/include/CesiumUtility/Math.h index c54ccb763..8a3dc49d1 100644 --- a/CesiumUtility/include/CesiumUtility/Math.h +++ b/CesiumUtility/include/CesiumUtility/Math.h @@ -3,6 +3,7 @@ #include #include +#include namespace CesiumUtility { @@ -233,6 +234,42 @@ class CESIUMUTILITY_API Math final { glm::vec(true); } + /** + * @brief Determines if two values are equal using an absolute or relative + * tolerance test. + * + * This is useful to avoid problems due to roundoff error when comparing + * floating-point values directly. The values are first compared using an + * absolute tolerance test. If that fails, a relative tolerance test is + * performed. Use this test if you are unsure of the magnitudes of left and + * right. + * + * @tparam T value value type. + * @tparam Q The GLM qualifier type. + * + * @param left The first value to compare. + * @param right The other value to compare. + * @param relativeEpsilon The maximum inclusive delta between `left` and + * `right` for the relative tolerance test. + * @param absoluteEpsilon The maximum inclusive delta between `left` and + * `right` for the absolute tolerance test. + * @returns `true` if the values are equal within the epsilon; otherwise, + * `false`. + */ + template + static constexpr bool equalsEpsilon( + const glm::mat<4, 4, T, Q>& left, + const glm::mat<4, 4, T, Q>& right, + double relativeEpsilon, + double absoluteEpsilon) noexcept { + for (glm::length_t i = 0; i < 4; ++i) { + if (!equalsEpsilon(left[i], right[i], relativeEpsilon, absoluteEpsilon)) { + return false; + } + } + return true; + } + /** * @brief Returns the sign of the value. *