fork download
copy 
  1. #include <algorithm>
  2. #include <cmath>
  3. #include <stdint.h>
  4.  
  5. #define PI_F 3.14159265358979f
  6.  
  7. template<typename T> class Vec2Base
  8. {
  9. public:
  10.  
  11.     Vec2Base() {}
  12.     Vec2Base(T a) : x(a), y(a) {}
  13.     Vec2Base(T a, T b) : x(a), y(b) {}
  14.  
  15.     T x, y;
  16. };
  17.  
  18. template<typename T> class Vec3Base
  19. {
  20. public:
  21.  
  22.     Vec3Base() {}
  23.     Vec3Base(T a) : x(a), y(a), z(a) {}
  24.     Vec3Base(T a, T b, T c) : x(a), y(b), z(c) {}
  25.  
  26.     friend Vec3Base<T> operator/(const Vec3Base& a, const Vec3Base& b)
  27.     {
  28.         Vec3Base<T> res;
  29.         res.x = a.x / b.x; res.y = a.y / b.y; res.z = a.z / b.z;
  30.         return res;
  31.     }
  32.  
  33.     friend T Dot(const Vec3Base& a, const Vec3Base& b)
  34.     {
  35.         return a.x * b.x + a.y * b.y + a.z * b.z;
  36.     }
  37.  
  38.     T x, y, z;
  39. };
  40.  
  41. typedef Vec2Base<float> Vec2f;
  42. typedef Vec3Base<float> Vec3f;
  43.  
  44. template<typename T>
  45. T SafeSqrt(const T& a)
  46. {
  47.     return std::sqrt(std::max<T>(0, a));
  48. }
  49.  
  50. template<typename T>
  51. T SignNum(const T& a)
  52. {
  53.     return std::signbit(a) ? T(-1) : T(1);
  54. }
  55.  
  56. Vec3f Normalize(const Vec3f& a)
  57. {
  58.     const float lenSqr = Dot(a, a);
  59.     const float len = std::sqrt(lenSqr);
  60.     return a / len;
  61. }
  62.  
  63. void CrashFooSample(
  64.     Vec2f       &aSlope,
  65.     const float  aThetaI,
  66.     const Vec2f &aSample)
  67. {
  68.     if (aThetaI < 0.0001f)
  69.     {
  70.         // Normal incidence - avoid division by zero later
  71.         const float sampleXClamped = std::min(aSample.x, 0.9999f);
  72.         const float radius = SafeSqrt(sampleXClamped / (1 - sampleXClamped));
  73.         const float phi = 2 * PI_F * aSample.y;
  74.         const float sinPhi = std::sin(phi);
  75.         const float cosPhi = std::cos(phi);
  76.         aSlope = Vec2f(radius * cosPhi, radius * sinPhi);
  77.     }
  78.     else
  79.     {
  80.         const float tanThetaI = std::tan(aThetaI);
  81.         const float tanThetaIInv = 1.0f / tanThetaI;
  82.         const float G1 =
  83.             2.0f / (1.0f + SafeSqrt(1.0f + 1.0f / (tanThetaIInv * tanThetaIInv)));
  84.  
  85.         // Sample x dimension (marginalized PDF - can be sampled directly via CDF^-1)
  86.         float A = 2.0f * aSample.x / G1 - 1.0f;
  87.         if (std::abs(A) == 1.0f)
  88.             A -= SignNum(A) * 1e-4f; // avoid division by zero later
  89.         const float B = tanThetaI;
  90.         const float tmpFract = 1.0f / (A * A - 1.0f);
  91.         const float D = SafeSqrt(B * B * tmpFract * tmpFract - (A * A - B * B) * tmpFract);
  92.         const float slopeX1 = B * tmpFract - D;
  93.         const float slopeX2 = B * tmpFract + D;
  94.         aSlope.x = (A < 0.0f || slopeX2 >(1.0f / tanThetaI)) ? slopeX1 : slopeX2;
  95.  
  96.         // Sample y dimension
  97.         // Using conditional PDF; however, CDF is not directly invertible, so we use rational fit of CDF^-1.
  98.         // We sample just one half-space - PDF is symmetrical in y dimension.
  99.         // We use improved fit from Mitsuba renderer rather than the original fit from the paper.
  100.         float ySign;
  101.         float yHalfSample;
  102.         if (aSample.y > 0.5f) // pick one positive/negative interval
  103.         {
  104.             ySign = 1.0f;
  105.             yHalfSample = 2.0f * (aSample.y - 0.5f);
  106.         }
  107.         else
  108.         {
  109.             ySign = -1.0f;
  110.             yHalfSample = 2.0f * (0.5f - aSample.y);
  111.         }
  112.         const float z =
  113.             (yHalfSample * (yHalfSample * (yHalfSample *
  114.             -(float)0.365728915865723 + (float)0.790235037209296) -
  115.             (float)0.424965825137544) + (float)0.000152998850436920)
  116.             /
  117.             (yHalfSample * (yHalfSample * (yHalfSample * (yHalfSample *
  118.             (float)0.169507819808272 - (float)0.397203533833404) -
  119.             (float)0.232500544458471) + (float)1) - (float)0.539825872510702);
  120.         aSlope.y = ySign * z * std::sqrt(1.0f + aSlope.x * aSlope.x);
  121.     }
  122. }
  123.  
  124. Vec3f CrashFoo(
  125.     const Vec3f &aVec3,
  126.     const float  aFloat,
  127.     const Vec2f &aVec2)
  128. {
  129.     const Vec3f vec3Stretch =
  130.         Normalize(Vec3f(aVec3.x, aVec3.x, std::max(aVec3.x, 0.0f)));
  131.  
  132.     float thetaWolStretch = 0.0f;
  133.     float phiWolStretch = 0.0f;
  134.     if (vec3Stretch.z < 0.999f)
  135.     {
  136.         thetaWolStretch = std::acos(vec3Stretch.z);
  137.         phiWolStretch = std::atan2(vec3Stretch.y, vec3Stretch.x);
  138.     }
  139.     const float sinPhi = std::sin(phiWolStretch);
  140.     const float cosPhi = std::cos(phiWolStretch);
  141.  
  142.     Vec2f slopeStretch;
  143.     CrashFooSample(slopeStretch, thetaWolStretch, aVec2);
  144.  
  145.     slopeStretch = Vec2f(
  146.         slopeStretch.x * cosPhi - slopeStretch.y * sinPhi,
  147.         slopeStretch.x * sinPhi + slopeStretch.y * cosPhi);
  148.  
  149.     Vec2f slope(
  150.         slopeStretch.x * aFloat,
  151.         slopeStretch.y * aFloat);
  152.  
  153.     return Vec3f(slope.x, slope.y, 1.0f);
  154. }
  155.  
  156. int32_t main(int32_t argc, const char *argv[])
  157. {
  158.     Vec3f vec3{ 0.00628005248f, -0.999814332f, 0.0182171166f };
  159.     Vec2f vec2{ 0.947231591f, 0.0522233732f };
  160.     float floatVal{ 0.010f };
  161.  
  162.     Vec3f vecResult = CrashFoo(vec3, floatVal, vec2);
  163.  
  164.     return (int32_t)vecResult.x;
  165. }
  166.  
Success #stdin #stdout 0s 3456KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/0I4OnD
language:
C++14 (gcc 8.3)
created:
9 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!