Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
22 commits
Select commit Hold shift + click to select a range
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions CHANGES.md
Original file line number Diff line number Diff line change
Expand Up @@ -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:

Expand Down
127 changes: 90 additions & 37 deletions Cesium3DTilesSelection/include/Cesium3DTilesSelection/ViewState.h
Original file line number Diff line number Diff line change
Expand Up @@ -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 {
Expand All @@ -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.
Expand All @@ -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,
Expand All @@ -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.
Expand All @@ -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.
Expand All @@ -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;
}

/**
Expand Down Expand Up @@ -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<CesiumGeospatial::Cartographic>& 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<CesiumGeospatial::Cartographic> _positionCartographic;

const CesiumGeometry::CullingVolume _cullingVolume;
const glm::dmat4 _viewMatrix;
const glm::dmat4 _projectionMatrix;
};

} // namespace Cesium3DTilesSelection
111 changes: 92 additions & 19 deletions Cesium3DTilesSelection/src/ViewState.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -5,17 +5,21 @@
#include <CesiumGeometry/CullingResult.h>
#include <CesiumGeometry/CullingVolume.h>
#include <CesiumGeometry/OrientedBoundingBox.h>
#include <CesiumGeometry/Transforms.h>
#include <CesiumGeospatial/BoundingRegion.h>
#include <CesiumGeospatial/BoundingRegionWithLooseFittingHeights.h>
#include <CesiumGeospatial/Cartographic.h>
#include <CesiumGeospatial/Ellipsoid.h>
#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/trigonometric.hpp>
#include <glm/ext/vector_double4.hpp>
#include <glm/geometric.hpp>

#include <cmath>
#include <limits>
#include <optional>
#include <variant>

Expand All @@ -39,34 +43,81 @@ 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,
const glm::dvec3& up,
const glm::dvec2& viewportSize,
double horizontalFieldOfView,
double verticalFieldOfView,
const std::optional<CesiumGeospatial::Cartographic>& 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<double>::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<double>::infinity()),
viewportSize,
ellipsoid) {}

namespace {
template <class T>
Expand Down Expand Up @@ -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;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function is called hundreds of times per frame, and it has gotten a lot more expensive in this PR. Fundamentally, I think it should still just be a linear function of geometricError / distance, right? Can we go back to computing something like the SSE denominator?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function is called hundreds of times per frame, and it has gotten a lot more expensive in this PR. Fundamentally, I think it should still just be a linear function of geometricError / distance, right? Can we go back to computing something like the SSE denominator?

The screen space error is not a function of distance for orthographic projections and is a bit messy for skewed perspective projections. In its current form this function doesn't need special logic for those cases; they are just a property of the projection matrix.

I am going to go out on a limb and assert that the performance of this function is dominated by the cost of floating point division. Before, this function had one division. At first glance the new version does 8, but only 2 of the results are actually used. Looking at the emitted assembly code with -O2, it looks like gcc does a good job of eliminating unused values, to the point that it only does 2 divisions and 10 multiplications. In other words, the dot products that would produce x and z values of the projection aren't done, and neither are the multiplications by 0.0 and 1.0.

It could be more clear to write the code so that it only does the operations that are needed and leave the more general math in comments, but I don't think that's a blocker?

I notice that the comment isn't quite correct and I will correct it.

Copy link
Member

@kring kring Apr 29, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can't see any evidence MSVC is optimizing this in a RelWithDebInfo build. The disassembly of this function looks like this:

double ViewState::computeScreenSpaceError(
    double geometricError,
    double distance) const noexcept {
00007FFFFADCBC96  mov         rdi,rcx  
  // Avoid divide by zero when viewer is inside the tile
  distance = glm::max(distance, 1e-7);
00007FFFFADCBC99  maxsd       xmm7,xmm2  
00007FFFFADCBC9D  movaps      xmmword ptr [rax-38h],xmm8  
  // 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);
00007FFFFADCBCA2  movups      xmmword ptr [rsp+20h],xmm0  
00007FFFFADCBCA7  lea         rcx,[rax-58h]  
00007FFFFADCBCAB  movaps      xmm8,xmm1  
00007FFFFADCBCAF  xorps       xmm7,xmmword ptr [__xmm@80000000000000008000000000000000 (07FFFFBA9B580h)]  
00007FFFFADCBCB6  movsd       mmword ptr [rsp+30h],xmm7  
00007FFFFADCBCBC  movsd       mmword ptr [rax-80h],xmm6  
00007FFFFADCBCC1  call        glm::operator*<double,0> (07FFFFA32C460h)  
00007FFFFADCBCC6  xorps       xmm0,xmm0  
  glm::dvec4 errorOffsetNdc =
      projMat * glm::dvec4(0.0, geometricError, -distance, 1.0);
00007FFFFADCBCC9  movsd       mmword ptr [rsp+28h],xmm8  
00007FFFFADCBCD0  lea         r8,[rsp+20h]  
00007FFFFADCBCD5  movsd       mmword ptr [rsp+20h],xmm0  
00007FFFFADCBCDB  lea         rdx,[rdi+1C8h]  
00007FFFFADCBCE2  movsd       mmword ptr [rsp+30h],xmm7  
00007FFFFADCBCE8  lea         rcx,[errorOffsetNdc]  
00007FFFFADCBCED  movsd       mmword ptr [rsp+38h],xmm6  
00007FFFFADCBCF3  call        glm::operator*<double,0> (07FFFFA32C460h)  
  errorOffsetNdc /= errorOffsetNdc.w;
00007FFFFADCBCF8  movsd       xmm0,mmword ptr [rsp+48h]  

  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;
}
00007FFFFADCBCFE  lea         r11,[rsp+0B0h]  
  errorOffsetNdc /= errorOffsetNdc.w;
00007FFFFADCBD06  divsd       xmm0,mmword ptr [rsp+58h]  

  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;
}
00007FFFFADCBD0C  mov         rbx,qword ptr [r11+10h]  
  centerNdc /= centerNdc.w;
00007FFFFADCBD10  movsd       xmm1,mmword ptr [rsp+68h]  
00007FFFFADCBD16  divsd       xmm1,mmword ptr [rsp+78h]  

  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;
}
00007FFFFADCBD1C  movaps      xmm6,xmmword ptr [r11-10h]  
00007FFFFADCBD21  movaps      xmm7,xmmword ptr [r11-20h]  
  errorOffsetNdc /= errorOffsetNdc.w;
00007FFFFADCBD26  subsd       xmm0,xmm1  

  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;
}
00007FFFFADCBD2A  movaps      xmm8,xmmword ptr [r11-30h]  
00007FFFFADCBD2F  xorps       xmm0,xmmword ptr [__xmm@80000000000000008000000000000000 (07FFFFBA9B580h)]  
00007FFFFADCBD36  mulsd       xmm0,mmword ptr [rdi+38h]  
00007FFFFADCBD3B  mulsd       xmm0,mmword ptr [__real@3fe0000000000000 (07FFFFBB4F748h)]  
00007FFFFADCBD43  mov         rsp,r11  
00007FFFFADCBD46  pop         rdi  
00007FFFFADCBD47  ret  

The important bit in my mind is those two calls to glm::operator*. If it's calling the mat4 * vec4 multiplication operator (twice) then that is a bunch of extra floating point operations that a) weren't done before, and b) shouldn't be necessary.

On the other hand, I can't measure any particular overhead from this, at least on my desktop on Windows. I'm still a little worried about low powered handheld and AR/VR devices, but it's probably true that there are bigger fish to fry. 🐟🍳

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I took screenshots from the same viewpoint in main and in this PR, and flipped between them, and I can't see any difference in LOD whatsoever. So no issues there.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It appears that MSVC isn't inliining the matrix multiply, which is unfortunate; if those calls are inlinined, then it's a pretty basic optimization to notice that results are unused and eliminate their computation. MSVC has managed to do that for the vector divide-by-scalar operations.

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
Loading
Loading