TessellationPredicates.cpp

#include "Utilities/TessellationPredicates.h"
 
#define _USE_MATH_DEFINES
#include <cmath>
#include <algorithm>
 
glm::vec2 Cogs::Core::localSphereRadiusInScreenSpace(const glm::mat4& /*localToClip*/,
                                                     const glm::mat4& clipToLocal,
                                                     const glm::vec4& clipPosition,
                                                     const float localRadius)
{
  glm::vec4 locPos4 = clipToLocal*clipPosition;
  glm::vec3 locPos3 = (1.f / locPos4.w)*glm::vec3(locPos4);
 
  glm::vec4 locPosDispX4 = clipToLocal*(clipPosition + glm::vec4(clipPosition.w, 0.f, 0.f, 0.f));
  glm::vec3 locPosDispX3 = (1.f / locPosDispX4.w)*glm::vec3(locPosDispX4);
 
  glm::vec4 locPosDispY4 = clipToLocal*(clipPosition + glm::vec4(0.f, clipPosition.w, 0.f, 0.f));
  glm::vec3 locPosDispY3 = (1.f / locPosDispY4.w)*glm::vec3(locPosDispY4);
 
  // radius in x and y pixels
  return localRadius*glm::vec2(1.f / glm::distance(locPos3, locPosDispX3),
                               1.f / glm::distance(locPos3, locPosDispY3));
#if 0
  float screenSpaceRadius = 0.f;
 
  glm::vec4 locPos4 = clipToLocal*screenPosition;
  glm::vec3 locPos3 = (1.f / locPos4.w)*glm::vec3(locPos4);
 
  {
    glm::vec4 xp4 = localToClip*glm::vec4(locPos3 + glm::vec3(localRadius, 0.f, 0.f), 1.f);
    glm::vec4 xm4 = localToClip*glm::vec4(locPos3 - glm::vec3(localRadius, 0.f, 0.f), 1.f);
    if ((xp4.z > -xp4.w) && (xm4.z > -xm4.w)) {
      glm::vec2 p = (0.5f / xp4.w)*viewportSize*glm::vec2(xp4.x, xp4.y);
      glm::vec2 m = (0.5f / xm4.w)*viewportSize*glm::vec2(xm4.x, xm4.y);
      screenSpaceRadius = glm::max(screenSpaceRadius, glm::distance(p, m));
    }
  }
  {
    glm::vec4 yp4 = localToClip*glm::vec4(locPos3 + glm::vec3(0.f, localRadius, 0.f), 1.f);
    glm::vec4 ym4 = localToClip*glm::vec4(locPos3 - glm::vec3(0.f, localRadius, 0.f), 1.f);
    if ((yp4.z > -yp4.w) && (ym4.z > -ym4.w)) {
      glm::vec2 p = (0.5f / yp4.w)*viewportSize*glm::vec2(yp4.x, yp4.y);
      glm::vec2 m = (0.5f / ym4.w)*viewportSize*glm::vec2(ym4.x, ym4.y);
      screenSpaceRadius = glm::max(screenSpaceRadius, glm::distance(p, m));
    }
  }
  {
    glm::vec4 zp4 = localToClip*glm::vec4(locPos3 + glm::vec3(0.f, 0.f, localRadius), 1.f);
    glm::vec4 zm4 = localToClip*glm::vec4(locPos3 - glm::vec3(0.f, 0.f, localRadius), 1.f);
    if ((zp4.z > -zp4.w) && (zm4.z > -zm4.w)) {
      glm::vec2 p = (0.5f / zp4.w)*viewportSize*glm::vec2(zp4.x, zp4.y);
      glm::vec2 m = (0.5f / zm4.w)*viewportSize*glm::vec2(zm4.x, zm4.y);
      screenSpaceRadius = glm::max(screenSpaceRadius, glm::distance(p, m));
    }
  }
  return 0.5f*screenSpaceRadius;
#endif
}
 
 
float Cogs::Core::sphereSectorGeometricDistancePredicate(const glm::vec2 screenSpaceRadii,
                                                          const glm::vec2 viewportSize,
                                                          const float distancePixelTolerance,
                                                          const float minSectors,
                                                          const float maxSectors)
{
  glm::vec2 pixelRadii = 0.5f*viewportSize*screenSpaceRadii;
  float radius = glm::max(pixelRadii.x, pixelRadii.y);
 
  //        __________
  //     _//\         \_
  //    / /c \          \
  //   | \__q \ r        |
  //  | /i  \__\          |
  //  |/_____phi\         |
  //  |a         b        |
  //  |                   |
  //   |                 |  
  //    \_             _/
  //      \___________/
  //
  // Let phi be the angle of the sector, and r is the radius of the circle, and
  // let c be the coord of the sector. We want to limit the distance between the
  // circular arc of the sector and the coord c.
  //
  // By bisecting the sector, we get q, which is orthogonal to c. If i is the
  // intersection between c and q, we get a triangle [a,b,i] where the corner
  // between [a,i] and [i,b] is pi/2, and thus the triangle is right. Hence,
  // the length of [i,b] is r cos(phi/2),
  //
  //           r - r cos (phi/2) <= tolerance
  //              1 - cos(phi/2) <= tolerance/r
  //    2 acos( 1 - tolerance/r) >= phi
  //
  // The number of sectors is n = 2pi/phi,
  //
  //    2 acos( 1 - tolerance/r) = 2pi/n
  //                           n = pi/acos( 1 - tolerance/r).
  // 
  static const float pi = float(M_PI);
  float a = 1.f - distancePixelTolerance / radius;
 
  // Make sure that we're inside acos' defined range
  a = std::max(a, std::cos(pi / minSectors));
  a = std::min(a, std::cos(pi / maxSectors));
 
  float sectors = float(M_PI) / glm::max(std::numeric_limits<float>::epsilon(), std::acos(a));
 
  return glm::clamp(sectors, minSectors, maxSectors);
}
 
float Cogs::Core::sphereSectorGeometricAreaPredicate(const glm::vec2 screenSpaceRadii,
                                                     const glm::vec2 viewportSize,
                                                     const float areaPixelTolerance,
                                                     const float minSectors,
                                                     const float maxSectors)
{
  glm::vec2 pixelRadii = 0.5f*viewportSize*screenSpaceRadii;
  float radius = glm::max(pixelRadii.x, pixelRadii.y);
 
  //        __________
  //     _//\         \_
  //    /C/  \          \
  //   | / T  \ r        |
  //  | /      \          |
  //  |/_____phi\         |
  //  |                   |
  //  |                   |
  //   |                 |  
  //    \_             _/
  //      \___________/
  //
  // Given an angle phi, the area of the triangle T is (1/2) r^2 sin(phi).
  // The area of the circular sector C is (pi r^2)/(phi/(2 pi)).
  // The approximation error for phi is area(C)-area(T),
  //
  //   (pi r^2)*(phi/(2 pi)) - (1/2) r^2 sin(phi) <= tolerance
  //             pi*(phi/(2 pi)) - (1/2) sin(phi) <= tolerance/r^2
  //                               phi - sin(phi) <= 2 tolerance/r^2
  //
  // The Maclaurin series for sin(x) is
  //
  //     sum_{n=0}^{\intfy} \frac{(-1)^n}{(2n+1)!}x^{2n+1} = x - x^3/3! + x^5/5! ...
  //
  // This is an alternating series (-1)^n with a_n = x^{2n+1}/(2n+1)!. If a_n+1 < a_n,
  // we can use Leibnitz alternating series theorem to cap the error of a truncated series.
  // Thus,
  //
  // a_{n+1} = x^{2(n+1)+1}/(2(n+1)+1)!
  //         = x^{2n+3}/(2n+3)!
  //         = (x^2/((2n+3)(2n+2))) x^{2n+1}/(2n+1)!
  //         = (x^2/((2n+3)(2n+2))) a_n
  //
  // and therefore if x^2 < (2n+3)(2n+2) the series converge. If we want truncate from
  // n=2, we get x^2 < (4+3)(4+2) = 7*6 = 36, and we have convergence for all sane values
  // of x. Further, the error is then bounded by x^{2*3+1}/(2*3+1)! = x^7/5040. In our
  // case, maximum sane phi represents three slices, i.e. (2/3)pi, which bounds the
  // error within 0.0004156.
  // 
  // Then, back to the inequality. We replace sin(phi) with the two first terms
  // of the Maclaurin series
  //
  // phi - phi + phi^3/3! + R <= 2 tolerance/r^2,   |R| < (((2/3)pi)^7)/(7!).
  //             phi^3 <= 6*(2 tolerance/ r^2 + R)
  //
  // phi = cbrt( 16 tolerance/ r^2) )
  //
 
  float phi = std::cbrt(16.f*areaPixelTolerance / (radius*radius));
 
  float sectors = 2.f*float(M_PI)/phi;
 
  return glm::clamp(sectors, minSectors, maxSectors);
}