Orthographic

Cesium3DTilesSelection/include/Cesium3DTilesSelection/ViewState.h

@@ -8,9 +8,11 @@
 #include <CesiumGeospatial/Ellipsoid.h>
 
 #include <glm/mat3x3.hpp>
+#include <glm/mat4x4.hpp>
 #include <glm/vec2.hpp>
 #include <glm/vec3.hpp>
 
+#include <optional>
 #include <vector>
 
 namespace Cesium3DTilesSelection {
@@ -25,7 +27,8 @@ 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.
    *
    * @param position The position of the eye point of the camera.
    * @param direction The view direction vector of the camera.
@@ -49,6 +52,25 @@ 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 The view's Vulkan style reversed Z projection
+   * matrix
+   * @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}.
+   */
+  static ViewState create(
+      const glm::dmat4& viewMatrix,
+      const glm::dmat4& projectionMatrix,
+      const glm::dvec2& viewportSize,
+      const CesiumGeospatial::Ellipsoid& ellipsoid CESIUM_DEFAULT_ELLIPSOID);
+
   /**
    * @brief Gets the position of the camera in Earth-centered, Earth-fixed
    * coordinates.
@@ -164,6 +186,13 @@ class CESIUM3DTILESSELECTION_API ViewState final {
       const std::optional<CesiumGeospatial::Cartographic>& positionCartographic,
       const CesiumGeospatial::Ellipsoid& ellipsoid);
 
+  ViewState(
+      const glm::dmat4& viewMatrix,
+      const glm::dmat4& _projectionMatrix,
+      const glm::dvec2& viewportSize,
+      const std::optional<CesiumGeospatial::Cartographic>& positionCartographic,
+      const CesiumGeospatial::Ellipsoid& ellipsoid);
+
   const glm::dvec3 _position;
   const glm::dvec3 _direction;
   const glm::dvec3 _up;
@@ -176,6 +205,7 @@ class CESIUM3DTILESSELECTION_API ViewState final {
   const std::optional<CesiumGeospatial::Cartographic> _positionCartographic;
 
   const CesiumGeometry::CullingVolume _cullingVolume;
+  std::optional<const glm::dmat4> _projectionMatrix;
 };
 
 } // namespace Cesium3DTilesSelection

Cesium3DTilesSelection/src/ViewState.cpp

@@ -12,8 +12,11 @@
 #include <CesiumGeospatial/S2CellBoundingVolume.h>
 
 #include <glm/common.hpp>
+#include <glm/ext/matrix_double4x4.hpp>
 #include <glm/ext/vector_double2.hpp>
 #include <glm/ext/vector_double3.hpp>
+#include <glm/ext/vector_double4.hpp>
+#include <glm/geometric.hpp>
 #include <glm/trigonometric.hpp>
 
 #include <optional>
@@ -43,6 +46,36 @@ namespace Cesium3DTilesSelection {
       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
+
+/* static */ ViewState ViewState::create(
+    const glm::dmat4& viewMatrix,
+    const glm::dmat4& projectionMatrix,
+    const glm::dvec2& viewportSize,
+    const CesiumGeospatial::Ellipsoid& ellipsoid) {
+  return ViewState(
+      viewMatrix,
+      projectionMatrix,
+      viewportSize,
+      ellipsoid.cartesianToCartographic(positionFromView(viewMatrix)),
+      ellipsoid);
+}
+
 ViewState::ViewState(
     const glm::dvec3& position,
     const glm::dvec3& direction,
@@ -68,6 +101,24 @@ ViewState::ViewState(
           horizontalFieldOfView,
           verticalFieldOfView)) {}
 
+ViewState::ViewState(
+    const glm::dmat4& viewMatrix,
+    const glm::dmat4& projectionMatrix,
+    const glm::dvec2& viewportSize,
+    const std::optional<CesiumGeospatial::Cartographic>& positionCartographic,
+    const CesiumGeospatial::Ellipsoid& ellipsoid)
+    : _position(positionFromView(viewMatrix)),
+      _direction(directionFromView(viewMatrix)),
+      _up(0.0),
+      _viewportSize(viewportSize),
+      _horizontalFieldOfView(0.0),
+      _verticalFieldOfView(0.0),
+      _ellipsoid(ellipsoid),
+      _sseDenominator(0.0),
+      _positionCartographic(positionCartographic),
+      _cullingVolume(createCullingVolume(projectionMatrix * viewMatrix)),
+      _projectionMatrix(projectionMatrix) {}
+
 namespace {
 template <class T>
 bool isBoundingVolumeVisible(
@@ -205,7 +256,28 @@ 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);
+  if (this->_projectionMatrix) {
+    // 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 is actually the slant distance, and the real
+    // transform contains a term for an offset due to a skewed projection which
+    // is ignored here.
+    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.
+    return -ndcError * this->_viewportSize.y / 2.0;
+  } else {
+    const double sseDenominator = this->_sseDenominator;
+    return (geometricError * this->_viewportSize.y) /
+           (distance * sseDenominator);
+  }
 }
 } // namespace Cesium3DTilesSelection

CesiumGeometry/include/CesiumGeometry/CullingVolume.h

@@ -2,6 +2,8 @@
 
 #include <CesiumGeometry/Plane.h>
 
+#include <glm/ext/matrix_double4x4.hpp>
+
 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 createOrthoCullingVolume(
+    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

CesiumGeometry/include/CesiumGeometry/Transforms.h

@@ -95,6 +95,83 @@ 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 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);
+
+  /**
+   * The projection matrix returned by the following functions uses the Vulkan
+   * 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")
+   */
+
+  /**
+   * @brief Compute a Vulkan-style perspective projection matrix with reversed
+   * Z.
+   *
+   * @param fovx horizontal field of view in radians
+   * @param fovy horizontal 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.
+   *
+   * @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.
+   *
+   * X maps from -1 to 1, Y maps from 1 to -1, Z maps from 1 to 0.
+   *
+   * @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